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ABSTRACT 


An analytical approach to the solution of the whipping response of a submerged submarine 
to the pulsation of a nearby explosion bubble is presented. Particular aspects of energy 
dissipation (damping) are addressed in regards to their effects on the overall loading forces 
and dynamic response of the submarine. Particular damping mechanisms discussed 
include hull damping (hull material and structural), internal damping (internal structures 
and sloshing liquids), and external damping (hydrodynamic damping including wave 
radiation and viscous fluid effects). 


A flexible analytic model is developed, using the fundamentals of vibration and hydro- 
mechanics. A finite-element beam model ts used to represent the flexural structure of the 
hull. Internal structures are modeled as separate dynamic systems, using both discrete and 
modal superposition approaches. Onboard liquids are modeled using a mechanical- 
equivalent system based upon a potential flow solution of the liquid free-surface. External 
hydrodynamic forces are modeled using a modified Morison formulation, with fluid 
velocities and accelerations calculated based upon a popular explosion bubble model. The 
equations of motion for the composite dynamic systems are solved ın the time-domain 
using a modified Newmark direct time integration scheme, with iterations at each time 
step accomplished using a modified Newton-Raphson method. 


The analytic model is implemented on a personal computer, and used to numerically 
investigate various damping mechanisms associated with dissipation of whipping energy. 
Additionally, the analytic model is compared to explosion-induced whipping tests 
conducted on the U.S. Navy test platform "Red Snapper". The analytic model compared 
well to the experimental data for early-time response (through the 2nd bubble period), but 
tended to over predict response for later times. 


Potential weaknesses with the analytic model are discussed in light of the comparisons 
with experimental data, and recommendations for future analytic work in the area are 
made. 


Thesis Supervisor: Professor Alan Brown, Department of Ocean Engineering 
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vector gradient operator 

volume 

matrix of mode shape vectors, velocity potential function 
Rayleigh coefficient, non-dimensional bubble radius 
Rayleigh coefficient, viscous-frequency parameter, coefficient for 
Newmark algorithm, control variable for free surface effect 
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strain, control variable for bubble migration 

coefficient for Newmark algorithm, adiabatic constant for explosive 
free surface elevation 

mode shape vector, velocity potential 

non-dimensional rate of change of bubble depth 

coefficient of friction 
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1. INTRODUCTION 


1.1 Underwater explosions and ship whipping. 

Designers of naval ships have long been concerned with the effects of underwater 
explosions. Historically, most attention has been paid to the effects of contact charges and 
charges of near to intermediate standoff. Protection methods to these effects have 
typically included heavy side-protective systems to minimize possibility of hull rupture and 
shock-mounting of vital machinery and equipment to minimize damage caused by the high 
frequency shock wave resulting from the explosion. 

While hull rupture and internal shock wave effects are still considered of major 
importance to ship survivability, increasing attention has been paid in recent years to the 
whipping of ship hulls caused by charges of near to intermediate standoff. Whipping can 
be defined as the transient, flexing motion of the whole hull of a ship in its vibrational 
modes corresponding to the lowest natural frequencies!. Such motion, if of sufficient 
magnitude, can cause hull girder buckling, tearing, or other loss of hull girder strength. 

Numerous examples of catastrophic or near-catastrophic whipping can be cited 
throughout naval history. As recently as 1991, the U.S. Navy Guided Missile Cruiser 
USS Princeton (CG-59) struck a floating contact mine ın the Persian Gulf. Although the 
local hull rupture caused by the mine explosion (near the bow) was not a threat to the 
survival of the ship, the subsequent whipping of the hull caused near-catastrophic hull 
girder damage near the stern of the ship. 

When an underwater explosive detonates, the solid explosive material suddenly 
reacts, leaving behind gaseous products at very high temperature and pressure. Almost 
immediately, the shock wave, produced due to the sudden discontinuity in pressure, 
propagates radially, approximately with the speed of sound in water. This shock wave 
carries with it approximately one-third of the energy released during the reaction”. It is 
this shock wave that is typically responsible for localized hull rupture and damage to 
internal equipment. 

Left behind as the shock wave propagates through the fluid medium is a cavity of 
gaseous reaction products at high pressure. This cavity subsequently expands (as a 
bubble) to relieve the difference in pressure, accelerating the surrounding fluid. The 
bubble continues to expand beyond the point of hydrostatic equilibrium (due to the inertia 


of the surrounding fluid) until a point of dynamic equilibrium is reached. The bubble then 


!Hicks. A.N.. "Explosion Induced Hull Whipping". Advances in Marine Structures, Admiralty Research 
Establishment. Dunfermline. Scotland, UK. 1986, p. 390. 
bid.. p. 392. 





reverses, continuing to contract until dynamic equilibrium is again reached, where it 
quickly rebounds and again begins to expand. This bubble expansion and contraction 
continues until the energy of the reaction is fully dissipated. As the bubble rebounds, it 
greatly accelerates the surrounding water, generating a substantial pressure pulse 
(commonly known as the bubble pulse). This bubble pulse can impart significant loads on 
structures in the vicinity. 

Thus, for certain underwater explosions, a ship can experience loading, not only 
from the radiating shock wave, but also from the quasi-periodic bubble pulse. Although 
the whipping motion of a ship hull 1s characterized by free vibration of the hull in the 
water, it is this quasi-periodic loading which can reinforce the free vibration and magnify 
the response of the hull to whipping, particularly when the period of the bubble pulse is in 


the vicinity of the natural periods of hull flexure. 


p The importance of damping in ship whipping. 

The first, most easily observed effect of damping on structural vibration is in the 
decrease in amplitude during free vibration. Figure 1-1 illustrates the effect of damping 
on the free vibration of a cantilever beam. If the cantilever with some damping is 
deformed and then released, it will oscillate with a regular period, but the amplitude of the 
oscillation will decrease with time. Without damping, the cantilever would continue to 
vibrate without the decrease in amplitude. In fact, the free vibrations of damped structures 
will normally be reduced at a rate that may be used as a measure of the amount of 
damping in the system. A well known measure of damping, known as the /ogarithmic 
decrement, 5, is related to the ratio of the th to the (7+ N)th cycle amplitudes by: 

ö= nl An 
N AN 
where A, is the amplitude of the wth cycle and A,,,, 1s the amplitude of the (7+ N)th 


cycle. Similarly to a cantilever, the amplitude of free vibration of a whipping ship will 


) (1-1) 





decrease with time due to the presence of damping. 

A second, and perhaps more important effect of damping is on the steady state 
amplitude attained when a structure is excited by a harmonically oscillating force. As 
discussed in any text on mechanical vibrations, the maximum magnitude of the response of 
any structure subjected to harmonic excitation (or periodic excitation that can be 
represented as a Fourier series of harmonic terms) is determined by the relation of the 


frequency of excitation to the natural frequencies of the structure, as well as the stiffness, 








Cantilever 


Figure 1-1. Effect of damping (with time) on the free vibration 
of a cantilever. 


damping and mass of the structure (and, of course, the magnitude of the excitation). As 
shown in figure 1-2 (illustrated for a single degree of freedom system subjected to 
harmonic excitation), the response is predominantly controlled by stiffness when the 
excitation frequency (o) is much less than the natural frequency (@,) of the structure 
(elastic forces dominate the equations of motion due to the small acceleration and inertia 
forces). When the excitation frequency is much greater than the natural frequency, the 
response is predominantly controlled by mass (inertia forces dominate the equations of 
motion). When the excitation frequency is close to the natural frequency (near 
resonance), however, the elastic and inertia forces approximately balance, and the major 
demand on the excitation forces is to overcome system damping, implying that the amount 
of system damping controls the dynamic response. In this region, when damping is low, 
the response magnitude is high, and when damping is high, the response amplitude is low. 
Thus, the response near resonance can be very much greater than the static response (o 
approaching zero), particularly if the damping is very low. 

As discussed previously, the bubble pulse can generate low frequency, quasi- 
periodic loading on the ship. For certain underwater explosive combinations (explosive 
charge weight and standoff), the frequency of this loading can be on the order of the 
natural frequencies of the ship's low frequency flexural modes. Thus, hull flexural 
vibrations, excited by the quasi-periodic bubble pulse loading, can approach a resonant 
response in the low frequency modes. It is clear then that damping can have a significant 


influence on the magnitude of the whipping response of the hull to the bubble pulse. 
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Figure 1-2. Effect of damping (with frequency) on the steady state 
vibration of a single degree of freedom system. 





2. DAMPING MECHANISMS IMPORTANT IN SHIP WHIPPING 


2.1 Introduction. 

In dynamics, the word damping generally refers to the dissipation of energy during 
the motion of a body or structure. This energy eventually is converted to heat by friction 
or generates waves in the surrounding medium. When a structure, such as a ship, vibrates 
or moves in on ocean environment, energy can be dissipated within the structure itself, 
through internal components, or externally into the surrounding fluid Specific 
mechanisms of damping, which can be expected to be important in ship and submarine 
whipping, are discussed here and then incorporated into a model for submarine whipping 


in chapter 3. 


2.2 Hull damping. 

The loss of vibrational energy within a complex structure can be accomplished 
through several possible mechanisms, depending on the frequency of the vibration and 
geometry of the structure. In general, low frequency energy (consistent with ship 
whipping) is lost within a structure through friction, which in turn generates heat which is 
conducted, convected and radiated away. This friction can occur among the material 
particles or crystals (within the structural elements) as they move relative to one another 
during material deformation. This damping mechanism is referred to as material or 
hysteresis damping. Friction can also occur at the interface between structural elements 
as they slide relative to one another (1.e. at joints). This damping mechanism is often 
referred to as structural or dry friction damping. 

As will be discussed in greater detail subsequently, it is useful to define the hull 
damping in terms of equivalent viscous damping. By doing this, the damping forces 
associated with material and structure can be represented as being proportional to the 
velocity. The resulting equations of motion for the vibration of the hull could then be 
written in the form 

Mx+Cx+Kx=f (2-1) 
where M is the mass matrix of the hull, C is the equivalent viscous damping matrix, K is 
the stiffness matrix, x is the displacement vector, f is the external load vector, and 
superposed dots refer to time derivatives (x is velocity and X is acceleration). Although 
material damping is more representative of the energy dissipation within the material, and 
structural damping more representative of losses in joints, it will be noted that an 


equivalent viscous damping is assumed for most rea/ structural dynamic analyses. The 





precise nature of the damping is less important than modeling the correct energy loss per 


cycle. 


2.2.1 Material (hysteresis) damping. 

When materials are deformed, energy is absorbed and dissipated by the material 
due to friction during internal reconstruction of the micro and/or macro structure as the 
material deforms. This reconstruction ranges from crystal lattice to molecular scale effects. 
This process, referred to as material damping or hysteresis damping, is discussed in great 
detail by Lazan? and Nashif*. 

All real materials dissipate energy during cyclic deformation. Regardless of the 
precise physical mechanisms involved, energy dissipation can be highly nonlinear, and 
therefore detailed analysis can be very difficult. One approach toward quantifying the 
internal damping behavior of real materials is through the Aysteresis /oop, figure 2-1. The 
loop plots strain vs. stress during cyclic deformation of a material volume and is obtained 
from experimental measurements. Typically, hysteresis loops representative of 
construction metals such as steel are extremely thin, deviating little from a single line 
(unless the metal is strained well into the plastic range). Thus, for elastic ship hull 
whipping of a steel hull, the damping provided by material hysteresis would be expected to 
be small compared to other damping mechanisms discussed subsequently. 

The area inside the hysteresis loop is the specific damping energy, D, and is equal 


to the energy dissipated per unit volume of material per cycle 


DZ $o-de= Tor) (2-2) 


where OG is stress and £ is strain. Equation (2-2) implies that the rate of energy dissipation 
(dD/dt) is proportional to the strain rate. This is equivalent to the hysteresis damping being 
proportional to the strain rate. 

The specific damping energy is very small for most conventional structural 
materials, such as steel, for typical elastic stress-strain levels. High-damping alloys and 
viscoelastic materials have higher specific damping energies, and therefore are much better 


at dissipating energy. 


3Lazan. B.J.. Damping of Materials and Members in Structural Mechanics, Pergamon Press. New York. 


1968. 
“Nashif. A.D.. Jones. D.I.G.. and Henderson. J.P.. Vibration Damping. John Wiley and Sons, New York, 


1985. 
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Figure 2-1. Typical hysteresis loops. 


As discussed in the previous section, it is useful to represent hysteresis damping in 
terms of an equivalent viscous damping. An equivalent viscous damping coefficient for 
hysteresis damping can be found by considering the harmonic motion of an equivalent 
spring-viscous damper system (figure 2-2). For this system, the force F(t) necessary to 
cause a displacement x(t) is given by 

F(t) = kx(1) + ex(t) (2-3) 
where k is the spring constant, c is the viscous damping constant, and x 1s the velocity 
(dx/dt). For harmonic motion of frequency o and amplitude X, the displacement can be 


written in the form 


xt) SINT OL) (2-4) 
Inserting equation (2-4) into equation (2-3) gives 
F(t) » kXsin(ot) 4 cXo cos(ot) (2-5a) 


F(t) kx tco J X^ - (Xsinoty (2-5b) 
Bear ca x ex (2-5c) 


When F(t) is plotted against x(t), equation (2-5c) represents a closed loop similar to a 
hysteresis loop (figure 2-3). The energy dissipated by the damper in one cycle of motion 


is found by substituting equation (2-5a) into equation (2-2): 
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(2-6) 
where AW is the energy loss in the damper in one cycle, which is the equivalent to the 


specific damping energy, D, for hysteretic damping. 
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Figure 2-2. Equivalent spring-viscous damper system. 





Figure 2-3. Equivalent spring-viscous damper system. 


For real hysteresis damping, however, it has been found experimentally that the 
specific damping energy, D, is independent of frequency, but approximately proportional 
to the square of the amplitude of the strain^. Thus, the equivalent viscous damping 


coefficient, c., Must be inversely proportional to the frequency, which can be represented 


by 


eq °’ 


eq 


Co = L (2-7) 
0 


Rao. S.. Mechanical Vibrations. Addison-Wellesley Publishing Co.. Reading. MA. 1990, p. 102. 


[5 





where h is called the hysteresis damping constant. As discussed subsequently, the idea 
that the equivalent viscous damping constant for material hysteresis damping is inversely 
proportional to frequency is made use of in the practical modeling of material damping in 


real structures. 


2.2.2 Structural (dry friction) damping. 

In addition to damping within the material of a structure, damping also occurs at 
the interface between structural elements as they slide relative to one another. This type 
of damping is usually referred to as structural damping, slip damping, Coulomb damping, 
or simply dry friction damping (it amounts to Coulomb or "dry" friction forces). This 
type of damping is considered important when structural elements are joined non-rigidly 
by riveting, clamping, or bolting, where a moderate amount of clamping pressure allows 
for some slip between elements. As for material damping, structural damping 1s discussed 
in great detail by Lazan®, and Nashif”. 

Because of the potential geometric complexities, quantitative evaluation of 
structural damping in complex structures 1s impractical even when numerical techniques 
such as finite element analysis are used. However, most structurally significant joints 
making up a ship hull are rigid/y we/ded rather than bolted or riveted. Thus, for elastic 
ship hull whipping, ıt would be expected that damping provided by structural dry friction 
forces would be very small compared to other damping forces discussed subsequently. It 
should be noted, however, that for certain ship designs (particularly many submarines and 
submersibles), internal decking ıs often connected to the main structure of the hull through 
various types of sliding deck joints or suspension mechanisms. Thus, in some ship 
applications, it Is necessary to account for some degree of dry friction damping, 
particularly when considering potentially large amplitude hull flexural motions. 

In general, Coulomb's Law of Dry Friction states that when two bodies are in 
contact, the tangential force required to produce sliding is proportional to the normal 
force acting in the plane of contact. Once sliding has begun, the force remains constant. 
The dry friction force is given by 

[UN (2-8) 
where N is the normal force and u is the coefficient of friction which is characteristic for 
the interface of the two surfaces. Thus, the structural dry friction damping force at the 


interface of two sliding surfaces is in a direction opposite to the direction of the relative 


azan. op. cit. 
TNashif. op. cit. 





velocities of the two surfaces but is independent of the displacement or velocity: it 
depends only on the normal force between the two surfaces. 

The fact that the dry friction forces remain constant as long as sliding motion exists 
between two surfaces (and are not dependent on displacement or velocity) leads to the 
conclusion that, in each successive cycle (or oscillation) of a structure, the amplitude of 
the motion is reduced by an amount proportional only to ihe friction forces resulting from 
the sliding within the structure. For a simple one dimensional case, the reduction in 
amplitude with each successive oscillation will remain constant as long as relative motion 


exists. 


2.2.3 Modeling material and structural damping in complex structures. 

If the material and structural damping properties of a structure could be 
quantitatively determined throughout the structure, a finite element procedure could be 
used to include the distributed damping properties (as element damping matrices) in the 
complete dynamic solution of the structural motions. In practice, however, direct 
evaluation of distributed structural damping properties in complex structures does not 
lend itself well to the numerical solutions for structural dynamics (although distributed 
damping has been included in some simple finite element applications*). Material and 
structural damping (as previously defined) are typically expressed in terms of damping 
ratios established from experiments, rather than by means of an explicit damping matrix. In 
fact, for many applications, there is no need to express the damping explicitly by means of 
a damping matrix, but rather only in terms of the modal damping ratios (G,) which can be 
applied through the principle of modal superposition. This concept is often called modal 
damping. An extension of modal damping in which an explicit damping matrix is defined 
through a linear combination of the mass and stiffness matrices is called Rayleigh 
damping. 

Modal damping and modal superposition. The method of dynamic analysis 
using modal superposition is discussed in numerous texts (Clough and Penzien?, Smith?, 
for example) and is discussed briefly here only for continuity. Fundamentally, the original 
set of N coupled equations of motion (for an N degree of freedom system) are 


transformed to a set of N independent normal (modal) coordinate equations. The 


®Kaliske. M.. Jagusch. J.. Gebbeken. N.. and Rothert. H.. "Damping models in finite element 
computations". Structural Dvnamics-Proceedings of the 2nd European Conference on Structural 


Dynamics (Vol. I), Rotterdam, 1995, pp. 585-591. 
?Clough. R.W. and Penzien. J.. Dvnamics of Structures (second edition), McGraw-Hill. New York. 1993. 
10Smith. JW.. Vibration of Structures - Applications in Civil Engineering Design. Chapman and Hall. 


London. 1988. 
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dynamic response can be obtained by solving separately for the response of each normal 
(modal) coordinate and then superposing these to obtain the response in the original 
coordinate system. 

The equations of motion for a multiple degree of freedom dynamic systems can be 
written in matrix form: 

MX - Cx - Kx - f (2-9) 

where, in general, matrices M, C, and K all have non-diagonal terms (i.e. the equations of 
motion are coupled). The modal equations of motion are obtained from equations (2-9) 


by applying a modal transformation, 


x= Dq (2-10) 
where © is the matrix whose columns are the mode shape vectors 
D=(0,0,...0, | (En 


(; 1s the vector representing the shape of the ith mode), and q is a vector of normal or 
modal coordinates. Pre-multiplying each term by the transform of the th mode shape 
vector, yields an equation of motion for each mode (uncoupled from other modes) of the 
form: 
pene gr Ta Os (2-12) 

where: 

M,=6,M6, =modal mass coefficient 

C. =6.Co, =modal viscous damping coefficient 

K, 2$, Kó, 2 modal stiffness coefficient 

Q. = 0,f =modal force (2-13) 
Equations (2-13) for modal mass and modal stiffness come about because of the 
orthogonality property of the normal mode shapes, which can be extended to modal 
viscous damping if the orthogonality condition is also assumed to apply to the damping 
matrix (as will be discussed subsequently). If equation (2-12) 1s divided by the modal 
mass coefficient, the modal equations of motion may be written in an alternate form: 





i, +26,0,4, +020, - a (2-14) 


n 


where o , is the (modal) natural frequency, and G, is defined as the modal viscous 


damping ratio 





E, ==» (2-15) 





The (modal) natural frequencies (@ , ) and mode shape vectors (, ) are found by solving 
the generalized eigenproblem for the undamped free vibration of the dynamic system: 
Ko, 2 0;M$, (2-16) 

As discussed previously, it is often more convenient (and physically reasonable) to 
define the damping of a complex multiple degree of freedom system using a viscous 
damping ratio for each mode, rather than to try to explicitly evaluate coefficients of the 
damping matrix C, because the modal damping ratios can be determined experimentally 
(or estimated with adequate precision in many cases). 

With the (uncoupled) modal equations of motion defined, the modal displacements 
can be found, and the total displacement can be found by summing the modal 


contributions using equation (2-10): 
N 
x= > 6.9. (2-17) 
n=l 


Thus, the total response of a system modeled using modal superposition can be thought of 
as a sum of the responses of N independent single degree of freedom systems. 

Rayleigh (proportional) damping. Under some circumstances, it is desirable or 
practicable to use modal superposition to uncouple the equations of motion in solving for 
the dynamic response of a system. In other circumstances, it may be desirable to use the 
basic concept of modal superposition, but develop an explicit damping matrıx (C) which 
may be used to solve the complete equations of motion in the time domain. This can be 
done by applying what has become known as Rayleigh damping or proportional damping 
(named after Lord Rayleigh, who first suggested its use). Rayleigh or proportional 
damping is discussed in varying detail in numerous sources (Clough and Penzien!!, 
Bathe!2, for example). 

Under Rayleigh or proportional damping, it is assumed that the damping matrix 
can be adequately represented by making it a linear superposition of the mass matnx (mass 
proportional damping), the stiffness matrix (stiffness proportional damping), or both. By 
applying a linear superposition, the orthogonality condition which applies to the mass and 
stiffness matrices (as discussed previously) would also apply to the damping matrix. Thus, 


the damping matrix could be written in the form: 
C=aM+ PK (2-18) 


11Clough and Penzien, op. cit. 
l2Bathe, K.J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, NJ, 
1982. 
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where a and ß are proportionality constants. The damping associated with equation (2-18) 
may be recognized by evaluating the generalized modal damping associated with it by 


combining it with equations (2-13) and (2-15), and solving for the modal damping ratio: 


C. =0,C0, =4 [aM+BK]ó, =aM, +BK_ (2-19) 
PAE E (2-20) 
20,M, 20, 2 


Thus, the damping associated with any particular mode can be thought of as a linear 
combination of a factor which is inversely proportional to frequency (the mass 
proportional term), and a factor which is directly proportional to frequency (the stiffness 
proportional term). An important point to note is that the Rayleigh damping model is 
consistent with the discussion of material (hysteresis) and structural damping presented 
earlier. Thus, Rayleigh damping can be thought of as providing a good model to account 
for material and structural damping (for elastic whipping, as discussed in the previous 
sections). 

The relationship between damping ratio and frequency for Rayleigh damping can 
be illustrated as shown in figure 2-4. The cases of purely mass proportional damping 
(B = 0) and purely stiffness proportional damping (a — 0) are shown separately from the 
combined or general case. 

As illustrated in figure 2-4, the two Rayleigh damping proportionality constants, a 
and p, can be evaluated by solving two simultaneous equations if the (modal) damping 
ratios, G, and G,, associated with two specific (modal) frequencies, o, and o ,, are 
known (i.e. from experiment). Substitution of the two damping ratios and frequencies into 
equation (2-20), combining the resulting equations into matrix form, and carrying out a 


matrix inversion gives the solution of the two damping proportionality constants: 


ed MEE oc 
Ic EE Ioco ES 


EA s a (uon) 
B ue oc E: 


It is apparent from figure 2-4 that the use of Rayleigh damping provides for very 
high damping ratios for modes of very high frequencies (o much greater than O ,). The 
end result is that the responses of very high frequency modes are effectively eliminated by 
the high damping imposed by the Rayleigh damping. In the case of elastic ship whipping, 


where only the lower few modes are of concern (as discussed previously), this limitation 





can be considered acceptable as long as the higher frequency used for evaluation of the 
Rayleigh proportionality constants (@ , ) is set among the higher modes expected to 


contribute significantly to the whipping response. 


Damping ratio 







Cambined 





Stiffneea propartianal 


Moss proportional | 


=< ee Frequency 


Figure 2-4. Relationship between damping ratio and frequency for Rayleigh damping. 


2.3 Internal damping. 

Within any complex structure, such as a ship, energy can be lost (or simply 
transferred) to "separate" dynamıc systems which are in some form "attached" to the main 
structure. For ships, two broad categories could be considered under this category. 
Discrete internal masses or structures such as machinery and equipment connected 
indirectly to the hull through decking, machinery foundations, and shock or sound mounts, 
can act as dynamic subsystems, redistributing and absorbing dynamic energy transferred 
from the hull. Additionally, onboard /iquids which have free surface can act as dynamic 
subsystems, redistributing and absorbing dynamic energy by sloshing. Because of the 
large number of discrete internal structures and liquid storage tanks found onboard most 
ships, it is important to consider the mechanisms through which energy can be transferred 


and absorbed into these internal damping devices. 


2.3.1 Damping through discrete internal structures. 
Heavy machinery, electronic equipment, and other discrete masses onboard a ship 
can be considered to have a significant potential for dissipating hull whipping energy. 
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Figure 2-5 illustrates notional discrete masses, connected to the hull via mounts, internal 
decking, foundations, etc. As will be discussed subsequently, the dynamic response of 
more complex internal structures could be modeled as easily as individual masses by 
applying the principle of modal superposition (under the assumption of linear response). 
For purposes of this analysis, the term "discrete internal structure" will be used to refer to 
any mass or discrete structure that may be modeled as a discrete internal dynamic system 
in terms of transfer and absorption of hull whipping energy. While it is certainly not 
desired that heavy machinery or electronic equipment absorb large amounts of energy, it 
must be assumed that their ultimate connection to the hull provides the potential for the 


transfer of energy from the whipping hull. 
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Figure 2-5. Discrete masses connected to the hull via 
mounts, internal decking, foundations, etc. 


Whether or not discrete internal structures connected to a hull provide damping to 
hull whipping depends upon the sizes of the masses, the dynamic properties of the 
mechanisms through which they are connected to the hull (decking, foundations, mounts, 
etc.), the dynamic properties of the hull, and the frequency of hull whipping (frequency of 
excitation). In section 1.2, 1t was stated that the controlling mechanism in the response of 
a system to excitation (stiffness, damping, or mass) ıs determined by the relation of the 
excitation frequency (or frequencies) to the natural frequency (or frequencies) of the 
system. The maximum amount of energy is transferred when the frequency of excitation is 
equal to a natural frequency (at resonance). In the case of discrete internal structures 
connected to a whipping hull, it is evident that maximum energy would be transferred 
when the frequency of hull whipping is equal to a natural frequency of the system 
composed of the hull, the discrete internal structure, and the connection mechanisms. An 
understanding of the dynamic effects of discrete internal structures connected to a 
whipping ship can be made by first considering the effects of a single mass connected to a 
(rigid) vibrating structure, where the structure and mass are free to vibrate in only one 
direction (a two degree of freedom system), and then, using the principle of modal 


superposition, extending this to more complex arrangements of multiple masses connected 
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to a vibrating structure, such as a whipping ship hull (i.e. multiple degree of freedom 
systems). 

The two degree of freedom system (dynamic vibration absorber). A single 
mass connected to a much larger vibrating structure can be used to dissipate or absorb 
vibrational energy of the larger structure. In machinery and structural dynamics, this 
concept has gained practical significance in the application of the dynamic vibration 
absorber. The concept of the dynamic vibration absorber is presented by Den Hartog'?, 
and their design and application is discussed in detail in numerous other sources 
(Korenev!?*, Hunt!”, for example). A small auxiliary mass is connected to a harmonically 
vibrating structure by a spring of specified stiffness (the location of the mass is selected to 
have the maximum effect on the significant resonant mode of the vibrating structure). In 
the case where damping is provided by the connecting mechanism (e.g. by a dashpot), the 
auxiliary mass 1s called a damped dynamic vibration absorber or tuned mass damper. 
Treating the main structure and vibration absorber as a two degree of freedom dynamic 
system (figure 2-6), Den Hartog showed that the amplitude of harmonic vibration of the 
main structure could be reduced to zero at a specified excitation frequency (or reduced 
over a range of frequencies). This could be accomplished by "tuning" the properties of 
the vibration absorber (adjusting the mass, stiffness of the spring, and damping in the 
dashpot). The downside in the application of dynamic vibration absorbers appears in the 
potentially large amplitude response (displacement and acceleration) of the absorber 
masses. In the case of heavy machinery and electronic equipment onboard a whipping 
ship, this could translate to high displacements and accelerations on the delicate 
equipment, and large forces on the mounts and decking. 

For the simple dynamic vibration absorber attached to a large vibrating structure, 
the equations of motion for this simple two degree of freedom system can be derived by 
considering dynamic equilibrium of each of the masses separately. Figure 2-7 illustrates 
each of the masses and the dynamic forces acting upon them. From equilibrium, the 
equations of motion may be written: 

MX+KX=F..+c(x-X)+k(x-X) (for the main mass) (2-23) 
mx +c(x-—X)+k(x-X)=0 (for the absorber mass) (2-24) 
where capital letters refer to the main structure and lower case letters refer to the 
absorber, as shown in figure 2-6. These equations of motion could be written in an 
alternate form by defining the force exerted on the main structure by the absorber: 
13Den Hartog, J.P., Mechanical Vibrations. (Jh edition), McGraw-Hill. New York. 1956. 


l14Korenev and Reznikov. Dynamic libration Absorbers, John Wiley and Sons, New York. 1993. 
15Hunt. J.B.. Dvnanuc Vibration Absorbers. Mechanical Engineering Publications. London. 1979. 
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HMM E RX) (2-25) 
which, of course is the negative of the force exerted on the absorber by the main structure. 
Thus, the equations of motion could be written: 


MX+CX+KX=F.,+F,. (for the main mass) (2-26) 
F corner = EX - X) +k(x- X)=-m& (for the absorber mass) (2-27) 
Equations (2-26) and (2-27) form a set of 2 simultaneous differential equations (coupled) 


which may be solved under various conditions for displacements, velocities and 
accelerations. 
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Figure 2-7. Dynamic forces on each mass of a dynamic vibration absorber. 


Multiple degree of freedom systems. For more complex arrangements of 
internal structure, where many degrees of freedom would be required to properly model 
the dynamic response, the displacements, accelerations, and velocities would be vectors, 
and a larger number of simultaneous differential equations would be required. For 


multiple degree of freedom main and internal structures, the equations of motion could be 





written just as equations (2-23) and (2-24), but in matrix form. Thus, there would be one 
simultaneous equation for each degree of freedom. 

An alternative approach for multiple degree of freedom systems (particularly 
where the internal structure is complex and a large number of degrees of freedom would 
be required to adequately represent the dynamics), would be to use the principle of modal 
superposition under the assumption of linear response (as discussed in section 2 2.3) to 
approximate the response of the internal structure. By including only those resonant 
response modes of the main structure and the internal structure which would contribute 
significantly to the overall response, the complexity of approximating the response of the 
main structure is greatly reduced. Figure 2-8 depicts a notional complex internal 


arrangement of discrete masses with many degrees of freedom. 
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Figure 2-8. Complex internal arrangement of discrete masses with 
many degrees of freedom. 


Using modal superposition, an equivalent mechanical model of the dynamic 
behavior of the discrete internal structure may be made as shown in figure 2-9. Each 
mode of resonant response of the discrete internal structure is represented by an 
equivalent modal mass, modal damping and modal stiffness (mass m, represents that 
portion of the discrete internal structure responding as a rigid body with the main 
structure). A variation of the dynamic equilibrium equation (2-24) for the dynamic 
vibration absorber can be written using the principle of modal superposition. The equation 
of motion for the nth mode of the internal structure can be written: 

m, X, +c, (š, -X)+k (x= X)=0 (2-28) 
where x, is the 7th modal displacement of the internal structure (equivalent to q, in 
section 2.2.3) and X is the displacement of the main structure (hull). 

The equations of motion for a system made up of an internal structure attached to 


a large vibrating structure (the hull) can be derived by considering dynamic equilibrium of 





the main structure and the superposed modal masses of the internal structure. The forces 
applied to the main structure by the internal structure can be obtained by modal 


superposition. Figure 2-10 illustrates the dynamic equilibrium for the main structure. 
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Figure 2-9. Equivalent mechanical model of the dynamic behavior of 
a discrete internal structure (using modal superposition). 
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Figure 2-10. Dynamic equilibrium for the main structure. 


From dynamic equilibrium, the equation of motion for the main structure can be 
written: 


[MX+CX+KX]=F,, +F,, (2-29) 





where F,, ıs an external force applied to the main structure, and F_, is defined as the force 
exerted on the main structure by the internal structure. F_, can be determined by 
considering dynamic equilibrium of the internal structure using modal superposition, as 


shown in figure 2-11. By dynamic equilibrium of the internal structure: 
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Figure 2-11. Dynamic equilibrium of the internal structure using modal superposition. 





It should be noted that the internal structure rigid body case (n = 0) can be taken out of 
equation (2-30) by noting: 
N N 


Y [s.s -X) e, s, -)]- -2m.s, - -m& - Yom, - -m, «M. X), (x, - X)] 


n=0 n=0 


(2-31) 
With this simplification, the force exerted on the main structure by the internal structure 


can be written: 


Bu A+ le (x, =X) + k, (x, — X)] (2-32) 


Equations (2-29), (2-28), and (2-32) form a set of N= / simultaneous equations (coupled) 
which can be solved under various conditions for displacements, velocities and 
accelerations. 

It should be noted that, for cases where more than one degree of freedom is 
important for the main structure (e.g. the main structure moves horizontally as well as 
vertically), the above method could be extended considering the modal response of the 
internal structure in each of the degrees of freedom of the main structure. Thus, the main 
structure displacement would be represented by a vector (X), the internal structure modal 


displacement would be represented by a vector (x, ), and the forces (F,, and F.,, ) would 


be represented by vectors. 





An application of damping through discrete internal structures will be included in 
the development of a model for submarine whipping in chapter 3. 

2.3.2 Damping through sloshing liquids. 

Liquid sloshing has been the focus of considerable attention in the space industry, 
due primarily to concerns about fuel-rocket interaction forces, and recognition that liquid 
sloshing resonance can result in significant dissipation of energy and can, therefore, be 
capitalized upon to suppress vibrations in large structures!®. Similarly, methods to 
account for the effects of sloshing liquids on the natural frequencies and damping of 
offshore platforms have been developed within the marine industry". The occurrence of 
resonance, effects of nonlinearities, and contribution of devices such as baffles have been 
the object of numerous other studies. It can be projected that liquids in tanks onboard a 
whipping ship could undergo sloshing resonance, resulting in the damping of whipping 
energy. Fundamentally, any onboard liquid with free surface has the potential for 
absorbing whipping energy through liquid sloshing (liquid without free surface acts only as 
a rigid mass). Figure 2-12 illustrates a notional liquid storage tank, restrained to move 
with the whipping hull, and the resulting sloshing (or standing waves) generated due to the 


existence of the free surface of the liquid in the tank. 
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Figure 2-12. Tank with sloshing liquid. 


The liquid in a moving tank can be expected to respond in a variety of ways 
dependent upon many variables (tank geometry, directions of tank motion, amplitudes and 


frequencies of tank motion, properties of the liquid, etc.). Translational or pitching 


lóWelt. F., and Modi. V.J.. "Vibration damping through liquid sloshing", Diagnostics, Vehicle Dynamics 
and Special Topics: 12th Biennial Conference on Mechanical Vibrations and Noise. ASME, 1989. 
\7Vandiver, J.K., and Mitome. S.. "Effect of liquid storage tanks on the dynamic response of offshore 
platforms". Dynamic Analysis of Offshore Structures: Recent Developments, CML Publications. 


Southampton. 1982. 





motions of a tank lead primarily to /areral sloshing (usually antisymmetric liquid mode 
shapes). Tank motions normal to the equilibrium free surface lead primarily to vertical 
sloshing (usually symmetric liquid mode shapes). Under conditions of abrupt or sudden 
acceleration, /iquid impact with tank overheads and opposing bulkheads may provide 
significant loading on tank boundaries. The high accelerations associated with a whipping 
ship may require consideration of any or all of these to account for the transfer of 
whipping energy. 

Analysis of liquid sloshing can be carried out by numerous means. Researchers in 
the field of liquid sloshing have used several principle methods to determine the important 
dynamic characteristics desired for various applications!®. In many applications, potential 
flow models, using linearized and nonlinear free surface conditions, have been developed 
and used to approximate the fluid motion in various container shapes (Bauer!?, Hutton°, 
Welt and Modi?!, for example). In order to properly account for energy transfer and/or 
dissipation, incorporation of liquid impact loading and viscous effects has been 
accomplished using equivalent mechanical models in conjunction with potential flow 
models (Dalzell??, Vandiver and Mitome??, Bauer?*, for example). 

Potential flow model. Fundamentally, lateral and vertical sloshing can be thought 
of as the motion of standing waves within a moving rigid container or tank. Under the 
assumptions of incompressible and irrotational flow, the motion of the liquid in a tank can 
be approximated well by solving the potential flow problem for the wave motion subject to 
appropriate boundary conditions (all variables are functions of time). The ve/ocity 
potential function ® (defined by v = V ® where v is the fluid velocity vector and V is the 
vector differential operator) represents solution of the differential equation (Laplace 
equation) 

V’®=0 (within the fluid domain) (2-33) 


subject to the following boundary conditions: 


18 Abramson. H.N. (ed.). The Dynamic Behavior of Liquids in Moving Containers, NASA SP-106. 1966. 
!9Bauer. H.F.. "Theory of the Fluid Oscillations in a Circular Ring Tank Partially Filled with Liquid". 
NASA Technical Note D-557, 1960. 

20Hutton, R.E., "An Investigation of Resonant. Nonlinear, Non-Planar Free Surface Oscillations of a 
Fluid", NASA Technical Note D-1870. 19653. 

?! Welt and Modi. op. cit. 

22Dalzell. J.F.. "Liquid Impact on Tank Bulkheads", The Dynamic Behavior of Liquids in Moving 
Containers, NASA SP-106, H.N. Abramson (ed.). 1966. pp. 353-372. 

23Vandiver and Mitome. op. cit. 

24Bauer, H.F.. "Nonlinear Mechanical Model for the Description of Propellant Sloshing". 4/-L4 Journal, 
Vol. 4. No. 9. 1966. pp. 1662-1668. 
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(at y ^ 1 (on the free surface)) (2-34b) 


where n is the normal vector of the tank wall, Vp is the velocity of the tank wall normal to 
its surface, g is the acceleration due to gravity, and n is the free surface elevation. The 
free surface boundary condition, equation (2-34b), represents the complete nonlinear free 
surface boundary condition?. The linearized free surface boundary condition (linearized 
about the mean free surface) would be written: 


^- d ^«p 
= *g——-0 (at y=0) (2-34c) 
m 








The fluid velocity vector (v) can be obtained from (by definition of the velocity potential): 
v= VO (2-35) 
The hydrodynamic pressure (P, ) exerted on the tank wall by the sloshing liquid can be 


obtained from the unsteady Bernoulli equations: 





Pp = zb: E ivo | (2-36) 
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where p is the liquid density. Finally, the hydrodynamic force exerted on the tank wall by 
the sloshing liquid can be obtained by integration of the hydrodynamic pressure over the 


area of the tank wall: 
F, = IE .n)dS (2-37) 
S 


where F is the hydrodynamic force vector (on the tank wall) and S is the surface of the 
tank wall (below the free surface). Thus, for a specified tank geometry and motion, the 
hydrodynamic forces exerted by the sloshing liquid on the tank can be approximated using 
a potential flow model. 

Equivalent mechanical model. As mentioned previously, in order to properly 
account for energy dissipation, viscous effects at the tank wall must be included. While it 
is difficult to include viscous effects in the potential flow solution, viscous damping may 


be included in an equivalent mechanical model for the sloshing. In an equivalent 


25Newman, J.N.. Marine Hydrodynamics. The MIT Press. Cambridge. MA. 1986. p. 247. 





mechanical model, each of the sloshing modes (or modes of standing wave motion) is 
represented by an equivalent spring-mass-damper as shown in figure 2-13. In typical 


practice, linear viscous damping is usually assumed for the (modal) dampers in the 


mechanical model?6, however, nonlinear mechanical modeling techniques do exist2” 
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Figure 2-13. Equivalent mechanical model of a tank with sloshing liquid 
using modal superposition. 


In an equivalent mechanical model, the size of the equivalent masses for each of 
the resonant sloshing modes (m, , m, ,...) may be determined by requiring that the 
(hydrodynamic) force exerted on the tank wall by the equivalent mass-spring system, for 
the case of zero damping, be equal to the force calculated from the potential flow theory 
for each resonant sloshing mode (m, is the portion of the liquid in the tank which acts as a 


rigid body). 





[PEE OST = (2-38) 





?6Silverman, S.. and Abramson. H.N.. "Damping of Liquid Motions and Lateral Sloshing", The Dynamic 
Behavior of Liquids in Moving Containers, NASA SP-106. H.N. Abramson (ed.). 1966. pp. 105-144. 
Bauer., op. cit. 





where f” is the hydrodynamic force exerted on the tank wall by the rth sloshing mode. 
The equivalent modal stiffnesses (k,,k,.,...) may be determined from the natural 


frequencies (@ ,,) calculated from the potential flow theory for each resonant sloshing 


Ik È 
a = a = = IA (2-39) 
n 


Damping coefficients for each resonant sloshing mode can be determined 


mode: 





experimentally for a given tank configuration via a logarithmic decrement method (section 
1.2). The modal damping ratio for the 7th mode (G,) can be determined from the log 
decrement (6) measured for each mode by 


= ERE x Moa (for small damping) (2-40) 


T n i 
where the log decrement is defined by equation (1-1). The modal damping ratio is the 
ratio of the modal damping coefficient to the critical modal damping coefficient (c/c.), 
where the critical modal damping coefficient represents that value of damping where the 
motion would not oscillate, but would be non-periodic. The modal damping coefficients 
(c,,c,,...) can be determined from 
c, = 26,m,0, (2-41) 

Referring to figure 2-13, the equation of motion for the ith equivalent mass can be 
written in terms of the displacement, velocity, and the acceleration of the equivalent mass 
as follows: 

m,X, -c,(X, - X)o k,(x, - X) 20 (2-42) 
where X is the displacement of the tank, x, is the displacement of the 7th sloshing mode 
(equivalent mass), and velocities and accelerations are represented by superposed dots. 
The total hydrodynamic force on the wall of the tank can be found by summing the forces 


contributed by each mode: 


F, --Y mk, --mX- Y m,&, (2-43) 
n=0 


n=l 
Similar to the discrete internal structure, each of the resonant modes of liquid sloshing 
can be applied to the main structure (the hull). 
As an example, consider the lateral motion of a rigid structure with an integral 


laterally sloshing tank. The equivalent mechanical model could be as shown in figure 





2-14. The equation of motion for each of the equivalent liquid masses (modal masses) is 
given by equation (2-42), and the total hydrodynamic force on the wall of the tank due to 
the liquid sloshing is given by equation (2-43). The equation of motion for the main 
structure can then be found by adding the total hydrodynamic force into the equation of 


motion for the structure without the liquid: 
| MX +CX+KX]=F., +F, (2-44) 
where F.., is an external force applied to the main structure. The hydrodynamic force (F, ) 


may be written (following the same principles of modal superposition used for the internal 


structure): 
F, --Y mx, 2-mX-« Ye, (x, - X) -k, (x, - X)] (2-45) 
n=0 nl 
It can also be noted that m,, the portion of the liquid which acts as a rigid body, can be 
found from: 
Denen em, (2-46) 
n=0 n=1 


where M, is the total mass of liquid in the tank. 
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Figure 2-14. Equivalent mechanical model for the lateral motion of a rigid 
structure with an integral laterally sloshing tank. 
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In terms of real applications, such as ın the whipping of a ship hull, equation (2-45) 
may be truncated because the equivalent masses m, and modal displacements x_ grow 
smaller with » (the rate at which they grow smaller depends upon the geometry of the 
tank). Ali terms should be included for which the natural frequencies of the tank (o ,) are 
less than or equal to the natural flexural frequency of the overall structure (less than the 
fundamental flexural/whipping natural frequency). If equation (2-45) is truncated to 
include N modes, then equations (2-44), (2-45), and (2-46) form a set of N- / 
simultaneous equations of motion for the total dynamic system consisting of the main 
structure and the sloshing liquid. 

It should be noted that, for cases where more than one degree of freedom is 
important for the main structure (e.g. the main structure moves horizontally as well as 
vertically), the above method could be extended considering the modal response of the 
sloshing liquid in each of the degrees of freedom of the main structure. Thus, the main 
structure displacement would be represented by a vector (X), the sloshing liquid modal 
displacement would be represented by a vector (x, ), and the forces (F, and F.) would be 
represented by vectors. 

An application of liquid sloshing will be included in the development of a model for 


submarine whipping in chapter 3. 


2.4 External (fluid) damping. 

The loss of vibrational energy of a whipping ship into the surrounding fluid 
medium can be accomplished through several general mechanisms. If the ship is 
sufficiently close to the free surface, or its motion is of sufficient velocity and/or 
frequency, surface and pressure waves of non-negligible energy can be generated. This 
form of fluid damping is often referred to as wave radiation damping. Additionally, 
whipping energy can be lost through mechanisms associated with the viscosity of the 
surrounding fluid. These mechanisms can be lumped under the general heading of viscous 
fluid damping, but include specific mechanisms associated with time-dependent viscous 
boundary layer shear forces ("skin friction" drag) as well as boundary layer separation 
("form" drag). 

In general, the dominant hydrodynamic forces on a structure are determined 
primarily by the characteristics of the fluid flow around the structure. For a whipping ship, 
the relative fluid flow around the hull would be generally oscillatory in nature. For 
oscillatory flow, the following non-dimensional parameters are often used to characterize 


the flow: 





U,D 





Peak Reynolds' number, Rn = 


V 

p 

Keulegan-Carpenter number, KC = NS 
R > 
Viscous-frequency parameter, B = En DE 
KC vT 


where U is the characteristic relative free stream velocity (U, is the amplitude of the 
characteristic relative free stream velocity or maximum relative velocity), D is the 
characteristic dimension of the structure (1.e. the diameter), v is the kinematic viscosity of 
the fluid, and T is the characteristic period of oscillation of the re/ative fluid flow. 
Additionally, the flow can be influenced significantly by the geometry of the structure, free 
surtace effects, sea floor effects, roughness on the surface of the structure, and the 
Orientation of the structure relative to the flow. It may also be shown that for 
harmonically oscillating flow around a fixed circular cylinder, KC = 27A / D, where A is 
the amplitude of the oscillation of the fluid particle displacement. Thus, KC may be 
thought of as representing an amplitude of motion between a structure and the fluid 
relative the characteristic dimension of the body. 

These parameters, in addition to the other influencing factors listed above, 
determine the characteristics of the flow, and thus the major components of the 
hydrodynamic force on the structure. As will be discussed in greater detail subsequently, 
these characteristic parameters (1.e. the characteristics of the flow) effect the occurrence of 
boundary layer separation, and thus a major component of viscous fluid forces (as will be 
discussed subsequently, boundary layer separation forces tend to dominate over boundary 
layer shear forces for bluff bodies). For fully submerged bodies, if KC ts /arge (1.e. U, 1s 
large and D is small), then boundary layer separation is likely to occur, and the major 
hydrodynamic forces are likely to be associated with viscous drag effects (the so-called 
"drag-dominated regime"). If both Rn avd KC are /arge (1.e. U, is large but D ıs xot 
small), then the added forces associated with the inertia of the surrounding fluid (called 
fluid inertia forces) must also be considered in addition to the drag forces (the so-called 
"Morison regime"). If KC is small (Rn may still be moderately large), then boundary 
layer separation is not likely to occur and the only viscous forces would be due to 
boundary layer shear. If, in addition to KC being small, the structure is on or near a free 
surface then radiation of surface waves tends to dominate the hydrodynamic damping (the 
so-called "diffraction/radiation" regime). It is important to note that, for complex 


structures such as ships, different portions of the structure may be subject to different 





types of hydrodynamic loading (depending upon the /oca/ characteristics of the structure 


and the relative fluid flow). 


2.4.1 Wave radiation damping. 

Fundamentally, waves are radiated whenever a structure moves in a fluid medium. 
As mentioned above, when KC is small (1.e. when D is large), and the structure is on or 
near a free surface, then forces associated with wave radiation are generally important. In 
typical practice, panel or boundary element methods are used to solve the boundary-value 
problem in analyzing hydrodynamic forces due to incident, diffracted and radiated 
pressures. Boundary element methods are based upon potential flow theory, which 
neglects the effects of viscosity. However, there has been some effort in recent years to 
include some of the effects of viscosity within the boundary element method (as will be 
discussed subsequently). 

Potential flow theory and its application to the hydrodynamics of structures is 
discussed in great detail in numerous sources (Newman??, Chakrabarti??, Faltinsen??, for 
example). Under potential flow theory (assuming inviscid, irrotational flow), the fluid 
velocity within the fluid domain may be represented as the gradient of a scalar potential 
function 

CD .c® co 
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(2-47) 


where v is the fluid velocity vector, P(x,y,z,t) is the scalar velocity potential function 
(written here in terms of Cartesian coordinates), and V is the vector differential operator. 
The velocity potential must satisfy Laplace's equation 


‚. do do Bo 
ET, (2-48) 
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CX Cy Qu. 
within the fluid domain, and satisfy boundary conditions specific to the geometry of the 
application. For the case of a body of general shape moving near a free surface in a fluid 
of finite depth (figure 2-15), the velocity potential must satisfy the following boundary 











conditions: 





Be = V.n (onthe surface of the body) (2-49a) 


cn 


28Newman. op. cil. 
29Chakrabarti. S.K.. Hydrodvnamics of Offshore Structures. Springer-Verlag. London, 1987. 
30Faltinsen. O.M.. Sea Loads on Ships and Offshore Structures. Cambridge University Press. UK. 1990. 
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(at y = n (on the free surface)) (2-49c) 
where n is the normal vector of the body surface (pointing from the body into the fluid 
domain), V is the velocity vector of the body, g is the acceleration due to gravity, n 1s the 
free surface elevation above the mean, and h is the mean depth of the water. The free 
surface boundary condition, equation (2-49c), represents the complete nonlinear free 
surface boundary condition?!. The linearized free surface boundary condition (linearized 


about the mean free surface) would be written: 
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Figure 2-15. Body of general shape moving in fluid with free surface. 


3lNewman. op. cit.. p. 247. 





Additionally, a radiation condition must be established (as a boundary condition limit at 
infinity) in order to ensure that radiated waves propagate away from the body. For 
solution of equation (2-48) in the time domain. initial conditions are also required: 
p — 0 

cp 

—->0 (st>-®) (2-50) 

cl 
Solution of Laplace's equation (2-48) for the velocity potential, ®, is obtained by the 
method of integral equations using Green's theorem which can be written in the general 


form: 


f(ov:G-Gv"oJav = [ (®VG - GV@)- nds CSi) 
V S 
where V is defined as the volume bound by surface S, where S includes the body surface, 
the free surface, the bottom surface, and the surface at "infinity" (figure 2-16). An 
appropriate Green function (G) can be obtained which satisfactorily represents the source 
potentials and their derivatives for various boundary condition combinations (Wehausen 
and Laitone??, Newman??, for example). 

Once the velocity potential function is obtained, the fluid velocity vector in the 
fluid domain can be found by definition of the potential function as written in equation 
(2-47). The hydrodynamic pressure exerted on the body by the fluid is obtained from the 


unsteady Bernoulli equations: 





P = -ol a + vol’) (2-52) 
Cl». 2 

where P is the hydrodynamic pressure (a function of position) and p is the fluid density. 

The hydrodynamic force exerted on the body (or the panel of interest) by the fluid is 

obtained by integrating the hydrodynamic pressure over the submerged portion of the 


surface: 
F= || (P-n)dS (2-53) 
S 


where F is the force vector on the surface or panel, S is the submerged surface of the body 


or panel, and n is the direction normal of the panel surface. 


32Wehausen. J.V.. and Laitone, E. V.. "Surface Waves". Encvclopedia of Physics (Vol. 9). Springer- 


Verlag. Berlin. 1960. 
33Newman, J.N., "Algorithms for the free surface Green function". Journal of Engineering Mathematics. 


19. 1985, pp. 57-67. 
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Figure 2-16. Surfaces of integration for Green's theorem. 


For slender bodies (where length 1s much greater than width), application of 
two-dimensional diffraction/radiation potential theory has proven accurate and useful, 
particularly at higher frequencies?*. Several two-dimensional approaches are available in 
the literature, the most common of which is strip theory??. In strip theory, the slender 
body is divided into many thin slices or strips, each of which has approximately constant 
cross-section, so that two-dimensional potential theory could be applied to determine the 
approximate hydrodynamic forces on each strip independently (figure 2-17). For rıgid 
body hydrodynamics applications, the hydrodynamic forces on each strip could then be 
integrated along the length of the body to determine the total hydrodynamic force on the 
body. Alternatively, similar to panel methods, the hydrodynamic forces on each "strip" 
can be resolved and integrated into the fluid-structure interaction dynamic problem as 
external loads. 

Thus, the motions of a /arge structure on or near a free surface (where viscosity 


effects can be ignored) could be obtained using the complete potential flow theory 


34Lewis, E.V. (ed.). Principles of Naval Architecture (Vol. III - Motions in Waves and Controllability), 
Second Revision, SNAME. Jersey City, NJ. 1989. 

35Salveson. N.. Tuck. E.O.. and Faltinsen. O.. "Ship Motions and Sea Loads". SVAME Transactions, Vol. 
78, 1970, pp. 250-279. 





(linearızed or fully nonlinear). It computes the forces exerted on the structure by the 
incident pressure (the external "loading" force - the bubble pulse for a whipping ship 
subjected to an underwater explosion bubble), the diffraction pressure (scattering of the 
incident pressure from the "rigid" structure), and the radiation pressure (wave radiation 
from the structure as it moves in each of its degrees of freedom). Thus, the total velocity 
potential (Dd) at any point in the fluid in the presence of a moving structure (either three- 
dimensional if panel methods are used, or two-dimensional if strip theory is used) can be 
written as a superposition of the velocity potentials due to incident, diffracted, and 


radiated pressures as follows: 


M MN 
P= O° +> +E (2-54) 
i 


iz] jel 
where ®° is the velocity potential due to the incident pressure, P% is the velocity potential 
due to the diffracted or scattered pressure from the /th panel (where M is the total number 
of panels making up the total surface S), and ®,, ıs the velocity potential due to the 


radiated pressure caused by motion of the ¡th panel in the jth direction (where N is the 


number of degrees of freedom for each panel) . 





Strip motion 


Figure 2-17. Two-dimensional strip theory. 
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For most marine applications, rigid body translations are ignored and the solution 
of the equations of motion is carried out in the frequency domain, assuming harmonic 
motion of the structure centered around an equilibrium position. The incident and 
diffracted velocity potentials are computed on the assumption that the structure is fixed in 
its equilibrium position, subjected only to the given external loading (usually wave 
loading). The radiated velocity potentials are then computed assuming that the structure 
undergoes one simple harmonic motion at a time (i.e. at each frequency in each degree of 
freedom) in calm water (see Chakrabart1**, Faltinsen*”, for example). For assumed simple 
harmonic motion of the form: 

ee (2-55) 
where x is the displacement vector, X is the displacement magnitude vector, and o is the 
frequency, the total hydrodynamic force on the structure (calculated with potential flow 
theory) has the form: 

F --o0'X,M,, *ioX,C,, (2-56) 
in which F, is the force on the structure in the ith degree of freedom, and X, is the 
displacement of the structure in the th degree of freedom. M _ ıs called the added mass 
matrix (hydrodynamic added mass in the wth degree of freedom due to unit displacement 
in the mth degree of freedom), and thus the first term on the right hand side represents the 
hydrodynamic forces which are in phase with acceleration. C... is called the radiation 
damping matrix (hydrodynamic damping in the wth degree of freedom due to unit 
displacement in the mth degree of freedom), and thus the second term on the right hand 
side represents the hydrodynamic forces which are iv phase with velocity. 

For applications where simple harmonic motion about a fixed equilibrium position 
IS of a reasonable assumption (e.g. when rigid body translations are important), the 
separate accounting for added mass and radiation damping forces as above is normally not 
possible (except in the case of steady rigid body translation). For cases with unsteady 
rigid body translation in addition to harmonic vibration (e.g. whipping of a submerged 
submarine subjected to an underwater explosion bubble), the calculation of the total 
hydrodynamic forces must be accomplished in the time-domain. Application of potential 
flow theory in the time-domain to predict genera! motions of structures is not new and 
fundamentally straight forward, but has had limited application for complex geometries 


and loadings because of limited computational resources. Discussions of direct time- 


36 Chakrabarti, op. cit. 
37Faltinsen. op. cit. 
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domain solution techniques are presented by various authors for various applications 
(Beck and Magee?®, Ogilvie*’, for example). 
2.4.2 Viscous fluid damping. 

As mentioned previously, dynamic energy can be lost into the surrounding fluid 
medium through forces associated with fluid viscosity. First, viscous shear forces occur 
within the viscous boundary layer surrounding the moving structure (so-called "skin 
friction" drag). Second, separation of the boundary layer from the structure leads to 
additional forces on the moving structure (so-called "form" drag). Together, boundary 
layer shear forces and boundary layer separation forces constitute the viscous fluid forces. 

The extent to which each of these forces 1s important is a function of the geometry 
of the body, and the relative fluid motion around the body (characteristic amplitude of 
relative fluid motion and characteristic frequency). Where geometry and fluid flow around 
the structure are such that the boundary layer remains attached (e.g. streamlined bodies at 
low angles of attack), boundary layer shear forces tend to dominate. However, for "bluff" 
bodies such as circular cylinders and "flat" plates, boundary layer separation readily occurs 
which tends to produce forces greatly exceeding the boundary layer shear forces. 
Numerous experiments have been conducted on various geometries giving indication 
under which conditions each of the forces would be significant, and estimation of the 
forces. 

The largest body of data exists for harmonic transverse oscillation of horizontal 
cylinders for which two dimensional flow around the body can be assumed. As discussed 
previously, flow around a two-dimensional body is typically characterized by several non- 
dimensional parameters, including the Keulegan-Carpenter number and the Reynolds' 
number (or alternatively the viscous-frequency parameter, B). Anaturk* and Anaturk et. 
al.*! present regions in the KC - B plane where two dimensional flows possess different 
characteristics. For example, when KC < | and B < 10° — 10° (Rn < 10° — 107), the 


boundary layer flow typically remains attached for the case of a smooth circular cylinder. 


38Beck. R.F., and Magee. A.R.. "Time-domain Analysis for Predicting Ship Motions". Dynamics of 
Marine Vehicles and Structures in Waves, Elsevier Science Publishers, B.V.. 1991, pp. 49-64. 

39Ogilvie. T.F.. "Recent Progress Toward the Understanding and Prediction of Ship Motions". 
Proceedings of the Fifth Symposium on Naval Hydrodynamics. Office of Naval Research. Washington, 
DC, 1964. pp. 3-128. 

40 Anaturk. A.R.. "An experimental investigation to measure hydrodynamic forces at small amplitudes and 
high frequencies". -tpplied Ocean Research, 13, 4, 1991, pp. 200-208. 

+1 Anaturk, A.R., Tromans. P.S., van Hazendonk. H.C.. Sluis. C.M.. and Otter. A., "Drag forces on 
cylinders oscillating at small amplitude: a new model". Journal of Offshore Alechanics and Arctic 
Engineering. 114. 1992, pp. 91-103. 
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As KC and ß (or Rn) increase, the turbulence in the boundary layer increases and flow 
separates. 

For the case where the two dimensional flow does not separate, the flow (and 
resulting forces) can be satisfactorily predicted using an attached laminar boundary layer 
theory (Newman", Stuart*, Batchelor**, for example). This theory predicts the drag 
forces resulting from the combined normal and tangential stresses in the fluid using 
Navier-Stokes equations. 

For a body of general shape, the forces associated with separated flow, and the 
point at which transition to separated flow occurs (i.e. the stability of the boundary layer) 
is dependent upon many factors. Many attempts have been made to solve separated flow 
around marine structures numerically (generally limited to two dimensional flow). 
Examples of various "state-of-the-art" methods include: 

(a) Single vortex method (Faltinsen and Sortland*>) 
(b) Vortex sheet model (Faltinsen and Pettersen?) 
(c) Discrete vortex method (Sarpkaya", Clements?*) 
(d) Vortex-in-cell method (Stansby and Dixon?) 
(e) Combination of Chorin's method and vortex-in-cell method (Chorin?9, 
Smith and Stansby?!) 
(f) Navier-Stokes solvers (Lecointe and Piquet??) 
Additionally, attempts have been made to include vortex shedding models in numerical 


solutions for potential flows and thus combine wave radiation with the effects of viscosity 


"Newman. op. cit. 

43Stuart. J.T.. "Double boundary layers in oscillatory viscous flow", Journal of Fluid Mechanics, 24. 
1966, pp. 673-687. 

#Batchelor. G.K.. An Introduction to Fluid Dynamics. Cambridge University Press. London. 1973. 
+>Faltinsen, O.M., and Sortland, B., "Slow-drift eddy making damping of a ship". Applied Ocean 
Research, 9. 1. 1987. pp. 37-46. 

*éFaltinsen, O.M., and Pettersen. B.. "Application of a vortex tracking method to separated flow around 
marine structures". Journal of Fluids and Structures, 1. 1987, pp. 217-237. 

+7Sarpkava. T., "Computational Methods With Vortices - The 1988 Freeman Scholar Lecture". Journal of 
Fluids Engineering, Vol. 111, March 1989, pp. 5-45. 

48Clements, R.R.. and Maull, D.J.. "The representation of sheets of vorticity by discrete vortices". 
Progress in Aerospace Science, Vol. 16. No. 2, 1975, pp. 129-146. 

49Stansby, P.K., and Dixon, A.G., "Simulation of flows around cylinders by a Lagrangian vortex scheme", 
Applied Ocean Research, Vol. 5, No. 3. 1985, pp. 167-178. 

°Chorin. A.J.. "Numerical study of slightly viscous flow". Journal of Fluid Mechanics, 57, 1973. 

pp. 785-796. 

31Smith. P.A.. and Stansby. P.K.. "Impulsively started flow around a circular cylinder by the vortex 
method", Journal of Fluid Mechanics, 194. 1988, pp. 45-77. 

32Lecointe, Y.. and Piquet. J. "Compact finite-difference methods for solving incompressible 
Navier-Stokes equations around oscillatory bodies". } on Karman Institute for Fluid Dynamics Lecture 
Series 1985-04, Computational Fluid Dynamics, 1985. 





(Yeung and Wu°?, Wong and Calısal°*, for example). While these methods have 
documented satisfactory agreement with experiment in some cases, they are generally not 
considered robust enough (numerical ınstabilities, etc.) to be applied with complete 
confidence to new problems where there is no guidance from experiments??. Thus, 
numerical methods to characterize separated flow are often used ¿n conjunction with 
experimental methods. 

Morison's equations. [n practice, when drag forces are important (and wave 
radiation may be ignored), and two-dimensional oscillatory flow can be assumed, a 
Morison's equation or modified Morison's equation can be used to predict hydrodynamic 
forces . Originally proposed by Morison et. al.°° for the prediction of wave forces on 
vertical piles which extend from the bottom through the free surface, the Morison 
equation has been modified to include application to cylinders of varying orientation under 
various loading conditions, including cases where the cylinders are free to move under the 
imposed fluid loading. 

Morison's equation for a stationary cylinder in an oscillating flow field includes 
the effect of fluid inertia forces and assumes drag forces to be a quadratic function of 


velocity: 


| 
F = pAC,ü  -pDC ,u[u| 2 


where F is the hydrodynamic force (per unit length of cylinder - sımilar to "strip theory", 
the total force on the cylinder can be found by integrating over the length of the cylinder), 
A is cross sectional area, D is diameter, u is fluid velocity (ù is fluid acceleration), 
coefficient C, is called the inertia coefficient (the first term being proportional to 


acceleration is equivalent to an inertia force), and coefficient C, is called the drag 


coefficient (the second term being proportional to velocity (squared) is equivalent to a 
drag force). Figure 2-18 illustrates the parameters used in Morison's equation. 

Morison's equation for an oscillating cylinder in a stationary fluid includes the 
effect of fluid added mass forces and again assumes drag forces to be a quadratic function 


of velocity, giving the reaction forces on the cylinder: 


55 Yeung. R.W., and Wu, C., "Viscosity effects on the radiation hydrodynamics of two-dimensional 
cylinders", Proceedings of the 10th International Conference on Offshore Mechanics and Arctic 
Engineering, ASME. 1991. pp. 309-316. 

5+Wong. L.H.. and Calisal. S.M.. "A numerical solution for potential flows including the effects of vortex 
shedding". Proceedings of the 1 1th International Conference on Offshore Mechanics and Arctic 
Engineering, ASME. 1992, pp. 363-568. 


55Faltinsen. op. cit. 
SéMorison. J.R.. O'Brien. M.P.. Johnson, J. W.. and Schaaf. S.A.. "The force exerted by surface waves on 


piles", Petroleum Transactions. AIME. 189, 1950, pp. 149-157. 
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"E ae 
E -pAC,X- -pDCix[s| (2-58) 


where x and x are the acceleration and velocity of the cylinder, C, is called the added 


mass coefficient, and C; is the drag coefficient for and oscillating cylinder in still water. 


E 
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Figure 2-18. Parameters used in Morison's equation 


Coefficients C_, C,, C,, and C, are normally experimentally determined for 
various geometries and loading conditions (usually as functions of KC and f or Rn), and 
have been determined in a large number of experimental studies (Bearman et. al.57, 
Sarpkaya??, Chakrabarti and Cotter??, for example). 

When a cylinder is free to move in a moving fluid, Morison's equation can be 


written in an extended form by combining equations (2-57) and (2-58): 


F - pAC, ü - Z pDC,u|u|- pAC,X - ZpDCA [s (2-59) 


This form is known as the dependent flow fields model*? and is obtained as the 
superposition of the two independent flow fields; a far field due to the water particle 


motion (i.e. the external loading) and relatively unaffected by the presence and motion of 


37Bearman. P.W.. Lin. X. W.. and Mackwood. P.R.. "Measurement and prediction of response of circular 
cylinders in oscillating flow". Proceedings of the 6th International Conference on the Behavior of 
Offshore Structures (BOSS), 1992. pp. 297-307, 

>8Sarpkava. T.. "Force on a cylinder in viscous oscillatory flow at low Keulegan-Carpenter numbers". 
Journal of Fluid Afechanics, Vol. 165. 1986, pp. 61-71. 

59Chakrabarti. S.K.. and Cotter. D.C.. "Hydrodynamic coefficients of mooring tower", Proceedings of the 
15th Annual Offshore Technology Conference, OTC 4496, 1983, pp. 449-458. 

60Chakrabarti. op. cit.. p. 187. 
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the structure, and a near field resulting from the motion of the structure. Coefficients C_ 
and C, may be obtained from experiments where the cylinder is held fixed in oscillating 
fluid, and coefficients C ,, and Cj may be obtained from experiments of an oscillating 
cylinder in otherwise still fluid. Alternatively, the coefficients could be determined from 


empirical relations (using KC and Rn, for example): 


E 





m 3) 


a i{ KC, =— Rn, = =p) (2-60) 





oT x.D 
(2-61) 


"e = (sc. = Rn = 


V 
where u, ,X, are the amplitudes of water particle and structure velocity, respectively, 
T, T' are periods of oscillation of water particle and structure, respectively, and subscripts 
Jand n refer to the far field and near field, respectively. 

An alternative to the independent flow fields model is the re/ative velocity model, 
in which forces are written in terms of the re/ative motion between the structure and the 


fluid, in which single coefficients are assumed to apply: 


F = pAC, (ü- X) x pAX ^ -pDC, (u - x)ju - | (2-62) 


In the relative velocity model, the coefficients may be determined using empirical relations 


(using KC and Rn, for example): 





T D 
C C= (c. - = - Rn Mm (2-63) 


V 
where subscript r refers to the relative motion, v, is the relative velocity between the 
structure and the fluid (v,, is the amplitude of the relative velocity), and T, is the 
combined period of the relative motion. The relative velocity model has seen extensive 
use in the analysis of hydrodynamic forces on drag dominated offshore structures (Laya et. 
al.°!, Spidsoe and Karunakaran®, for example). 
Thus, the Morison equation can be used under various conditions to predict the 


hydrodynamic forces on structures (or portions of structures) where drag forces and fluid 


6lLaya, E.J.. Conner, M.. and Sunder. S.S.. "Hydrodynamic forces on flexible offshore structures". 
Journal of Engineering Mechanics, ASCE. March 1984. pp. 433-448. 

62Spidsoe. N., and Karunkaran, D.. "Nonlinear effects of damping to dynamic amplification factors for 
drag-dominated offshore platforms". Proceedings of the 11th International Conference on Offshore 
Mechanics and Arctic Engineering, ASME, 1992, pp. 231-239. 
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inertia forces are important. An application of Morison's equation to the development of a 


model for submarine whipping will be made in chapter 3. 
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3. MODEL FOR SUBMARINE WHIPPING INCLUDING DAMPING EFFECTS 


3.1 Introduction. 

In order to gain some quantitative insight into the re/ative amount of damping 
which could be provided by each of the mechanisms discussed in the previous chapter, it is 
first necessary to develop a genera! model for submarine whipping, which may include arıy 
or all of the specific mechanisms. One way of doing this would be to develop a simplified 
model of each of the mechanisms based upon the fundamental principles of mechanics, 
some of which were used in the discussion in the previous chapter, and utilize a flexible 
formulation for the equations of motion which permits inclusion of se/ected mechanisms. 
Using this approach, any individual mechanism could be evaluated separately for its 


expected effect on total damping. 


The equations of motion. The dynamic behavior of any structure is described by 
the equations of motion, which are simply extensions of Newton's second law of motion. 


Fundamentally, an equation of motion can be written as a sum of forces which resist 


motion of the structure (i.e. the inertial forces (F „~ ), the damping or dissipative forces 


inertia 


(Finis; ), and the elastic or stiffness forces (F,,,..,.)), and forces which induce motion (1.e. 


the external or loading forces (F )). Thus, an equation of motion may be written in the 


loading 
general form: 


F +F +F =F 


damping loading (3 | ) 


Inertial forces are those forces through which kinetic energy is stored, and are 


inertial elastic 


dependent upon acceleration and a mass coefficient (F,...,,(M,x)). Elastic forces are 
those forces through which potential energy is stored and are generally dependent upon 
displacement and some elastic or stiffness coefficient (F.,,,.,.(K,x)). As discussed 
previously, damping forces are those forces related to the dissipation of energy, and 
depend upon the mechanisms through which energy is dissipated, generally dependent 
upon velocity and some damping coefficient (Fy.noine (C.X)). It should be noted that, in 
general, these dependencies are neither constant nor linear. 

The external loading forces include all forces which are "externally" applied to the 
structure. These forces may themselves be thought of as being inertial, dissipative, or 
elastic in nature. Thus, equation (3-1) may be rewritten by redefining the inertial, 
damping, and elastic terms such that they also include the "external" loading forces. For 
example, for a system consisting of a structure moving in a fluid, equation (3-1) may be 


rewritten 
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JS +F...=0 


| 
F damping elastic ( D 2) 


ınertial 


| . . . . . . 
where Faena represents a// inertial forces (inertia of the structure and inertia of the 
surrounding fluid), F eus represents a// damping or dissipative forces (energy dissipation 


within the structure and energy lost into the surrounding fluid), and F viene fepresents all 


lastıc 
elastic forces (structural stiffness and hydrodynamic buoyancy forces). 

In general, the terms associated with inertial forces, damping forces, and elastic 
forces may be written in any order as long as all forces are accounted for (and as long as a 
proper sign convention is maintained). In keeping with the classification of damping 
forces discussed in chapter 2, the equations of motion for a whipping ship might therefore 


be written in the form: 


Fun + E^ emal t P = Ö (3-3) 
where: 
Faun = Lena Al ne + De po. 
SA emal = Lens + D a en n 22] 
pm x Fes as Den D E hol (3-4) 
where F «ema would include a// external forces on the hull. Thus, the dynamics of a 


extema 


whipping ship may be thought of as a (coupled) superposition of the dynamics of the hull, 
any internal components, and the external fluid medium. Equations (3-3) and (3-4) 
represent one form of the complete, coupled, nonlinear equations of motion for a 
whipping ship. 

Using this approach to the equations of motion, each term may be represented by 
an independent formulation (which may involve specific assumptions of linearity for thar 
term), so long as the terms remain coupled via a single solution technique. As an example, 
the dynamic behavior of a slender ship /n/// (such as a modern submarine) might be 
calculated based upon a linear simple beam finite element model, while the external 
(hydrodynamic) forces might be calculated based upon a nonlinear three-dimensional panel 
method. As will be discussed subsequently, even if the hull itself can be assumed to 
respond linearly, the nonlinearities associated with the external (hydrodynamic) forces 
require the solution of the equations of motion in the time-domain. It is the existence of 
the coupling of the terms within the equations of motion which require the solution of the 


combined equation by a single numerical algorithm. 
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9 Formulation of the dynamic analysis procedure. 

In general structural dynamics, the selection of the most appropriate method to 
solve the equations of motion depends upon the type of loading on the structure, the 
extent to which the equations of motion can be uncoupled and linearized, and the 
existence of rigid body motions. Nonlinearities which must be addressed include 
geometric nonlinearities (large rotations), material nonlinearities (nonlinear material stress- 
strain relations), nonlinearities associated with the external loading (hydrodynamic 
nonlinearities), nonlinear boundary conditions, and so on. The extent to which 
nonlinearities can be "linearized" depends upon the assumptions made in developing the 
physical model for the motions of the particular structure. For ship whipping, even if 
material and geometric nonlinearities can be effectively linearized (by assuming linear 
elastic stress-strain relations and small rotations - reasonable assumptions for the case of 
elastic whipping of a stee/ ship hull), hydrodynamic forces can be expected to be highly 
non-linear (as discussed in chapter 2), and thus complete linearization of the equations of 
motion is generally or possible Additionally, for the general dynamic response of a 
whipping ship subjected to an underwater explosion bubble, it could be expected that there 
would be a significant potential for rigid body motion. For these reasons, the equations of 
motion must be solved directly in the time-domain (utilizing a time-stepping numerical 
algorithm). 

Implementation of time-domain analysis procedures for the dynamic response of a 
structure in a fluid medium is discussed in great detail in a large number of texts and 
journals articles. Bathe®?, Clough and Penzien**, and Argyris and Mlejnek°° provide 
particularly thorough discussions of the solution of nonlinear equations of motion by direct 
time integration methods. Direct time integration means that the equations of motion 
would be solved using a numerical step-by-step procedure. The term "direct" generally 
implies that there is no transformation of the equations into a different form prior to 
numerical integration, although using the formulation of the equations of motion discussed 
above, it is feasible to utilize a coordinate transformation, so long as it is conducted at 
each time step (this will be discussed in greater detail later). In essence, instead of trying 
to satisfy equilibrium at only one distinct time / (as for static equilibrium), direct time 
integration aims to satisfy dynamic equilibrium at a// discrete time intervals Av apart. This 


implies that, basically, dynamic equilibrium (which includes the effects of inertia, damping 


63Bathe. op. cit. 
63Clough and Penzien, op. cit. 
65 Argyris. J.. and Mlejnek. J.P.. Dvnanmics of Structures, Elsevier Science Publishers. The Netherlands. 
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and elastıc forces) is sought at discrete time points within the interval of solution. 
Typically, for direct time integration algorithms, the solution at the "next" required time, 
t* At, 1s calculated based upon the solutions at all previous times, and assuming an 
incremental variation in displacements, velocities, and accelerations. It is the form of the 
assumption on the variation of displacements, velocities, and accelerations within each 
time interval that determines the accuracy, stability, and required time of the solution 
procedure. 

For direct time integration methods, the incremental equation of motion can be 
written as an extension of equation (3-2) by subtracting the equilibrium equation at time t 


from that at time t+At: 


T = F. 


inertial inertial 


ee SE Ö (3-5) 


elastıc 


ESI EE 


damping damping 
As discussed above, inertial forces are dependent upon acceleration and a mass coefficient 
(F, 4, (M. X)), damping forces are dependent upon velocity and some damping coefficient 
( Fio; (C. X)). and elastic forces are dependent upon displacement and some elastic or 


stiffness coefficient (F.,,...(K,x)). Thus, if coefficients M, C, and K can be related to 


lastic 
acceleration, velocity and displacement (respectively), then equation (3-5) contains only 
three unknowns at time t+At (the solution at time t would be known) and may be solved if 
two of the three are known as relations to the third. 

The Newmark time integration method. The Newmark integration method® , 
which is one of the most reliable direct time integration schemes, and whose stability 
condition is well known, can be employed quite appropriately for this application (see also 
Bathe®’, Argyris and Mlejnek®’, and Doyle9?). In the Newmark method (often referred to 
as the Newmark B Method), iterative calculation is carried out at each time step and is 
moved to the next time step after an equilibrium condition is reached. 

The basic equations of the Newmark method are derived by assuming Taylor 
expansions for displacements (x) and velocities (x ): 


> 


ee 


: .. At EE r > 
-X FX ALEX, RX yo XX, JA (3-6) 


Kr 


x = x, E X,At DE X, )At (850) 


t+ At 


66Newmark. N.M.. "A method of computation for structural dynamics." American Society of Civil 
Engineers Transactions. Vol. 127, 1962. pp. 601-630. 
67Bathe. op. cit. 


68 Argyris and Mlejnck. op. cit.. pp. 485-488. 
69Doyle, J.F., Static and Dynamic Analvsis of Structures. Kluwer Academic Publishers. Boston. 1991. pp. 
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where ß and y are coefficients to evaluate contributions of remainders. Thus, the 
displacements and velocities are approximated based upon “assumed” accelerations at 
each time interval (thus requiring an iterative approach), and based upon the solution at 
the previous time interval. Newmark showed that y and f) are not entirely arbitrary and 
their selection determines the stability and accuracy of the numerical solution. Newmark 
found that unless y is greater than or equal to 0.5, spurious damping is introduced into the 
equations of motion. Additionally, unless B is greater than or equal to 0.25, there is an 
incorrect maximum velocity response. Further, with y = 0.5 and B = 0.25, the scheme 
leads to unconditionally stable solutions. 

To implement the Newmark method in the analysis of the dynamics of a structure 
subject to nonlinear loading forces (such as a whipping submarine), the equations of 
motion are first put in the form: 

EN. -F loading 


where F refers to all those forces associated with the structure itself which resist 


resisting 


motion (i.e. F,,, as defined in section 3.1) and F,.,,,,. refers to all (other) forces which 


as defined in section 3.1). Defining F for a 


resisung 


+F 


external 


induce motion (1.e. F 


int ernal 
structure as: 
Fo. =Mx+Cx+Kx (3-9) 


resisting 
the equations of motion for a structure can be written in the general form: 


Mx +Cx+Kx =F (3-10) 


loading 
or in incremental form: 

MG. -X)RCO TX) FKG S 7X)T Eun Bas, — 6-11) 

From equation (3-10), the acceleration required for dynamic equilibrium to exist at 


time t+At can be written: 


ae = M" Lu ioe E eN 2 Res ) (3-12) 


Because F is generally also dependent upon the accelerations (and displacements and 


loading 


velocities by equations (3-6) and (3-7)), it is evident that an iterative procedure must be 


implemented at each time t+At in order to obtain equilibrium by equation (3-12). 





Newmark’? proposed a general method for the solution of structural dynamics 
which may be summarized as follows: 
IE Assume initial values of displacement, velocity, and acceleration (i e. at 


each node) for time t = 0. 


ID 


For each time step: 

(a) Assume values of acceleration for the current time step. 

(b) Compute displacement and velocity for the current acceleration 
(from equations (3-6) and (3-7). 

(c) For the current displacement, velocity, and acceleration, compute 
the loading and resisting forces in accordance with the physical 
model for the particular application. 

(d) Compute the acceleration required to obtain dynamic equilibrium 
with these loading and resisting forces (from equation (3-11)). 

(e) Compare the required acceleration with the assumed acceleration 
for the current time step. If these are the same (1.e. within 
tolerance), the calculation for this time step is complete. If these 
are different, repeat the calculation (1.e. perform another iteration 
through steps (a) - (e)) with a different value of assumed 
acceleration (the derived acceleration could be used as the new 
assumed acceleration). 

Move onto the next time step (t-At), assuming the final acceleration for the 


u 


last time step (t). 

It should be briefly noted here that the modeling of internal damping devices such 
as structures or liquids requires only the addition of N equations of motion to the set 
defined above (for an internal damping device with N participating modes of response). 
Thus, in addition to the requirement for equations (3-12) to be satisfied at each time t+At, 


the following equations must also be satisfied at each time t+At: 


X == 2 D (3-13) 
m 


n 


where subscript n is for each mode of the ınternal damping device (up to mode N). Thus, 
if there are M degrees of freedom required to represent the displacements of the hull, and 
N (modal) degrees of freedom required to represent the (modal) displacements of the 
internal device, then a total of M - N simultaneous equations must be satisfied using the 


form of equations (3-12) and (3-13). 


70Newmark. op cii CDD: 606-611. 





To implement the Newmark method in the analysis of the dynamics of a structure 
subject to loading that is highly nonlinear (such as a whipping submarine subject to 
nonlinear hydrodynamic forces), a slight modification to the above method must be made. 
Because of the highly nonlinear loading forces, the equations of motion must be solved at 
each time step using an iterative process (such as a Newton-Raphson process) which 
"guarantees" convergence. Without such an iterative process, the solution may not 
converge and equilibrium may not be (numerically) obtained. 

Suzuki and Yoshida?! discuss application of the Newmark method with solution at 
each time step by an iterative Newton-Raphson method. In the case of elastic ship 
whipping with nonlinear hydrodynamic loading, the specific formulation of the Newton 
iteration may proceed by using equations (3-6), (3-7) and (3-11), and expressing 


acceleration and velocity at time t+At in terms of displacement at time t+At 


n ] l l 

Xa = — (Xy -x,)- Xi le: (3-14) 
sc mE Dee END 28 J^ 

. m UE Yale a 

X e = uch XA 3-15 


In order for equilibrium to exist at time t+At, equations (3-11), (3-14) and (3-15) must be 
simultaneously satisfied. Residuals (which must vanish for equilibrium to exist) may be 
defined for equations (3-11), (3-14) and (3-15): 


i x m esr x IA ER + |] F ES 
eM sse x) um | 5| das. E Lk 1 sat] peer thee 
(3-16) 
> : Var le : 
P run mE u x,) (i i | EJ 3-17) 
. | d | i 
pal ice ae t+At -x)- s X [x (3-18) 


The Newton iteration is carried out such that the unsatisfied equilibrium at iteration / is to 
be satisfied at iteration /+/. In other words, all residuals at iteration / - / must 


simultaneously vanish. To obtain expressions for residuals at iteration step /~ /, the 


7lSuzuki, H.. and Yoshida. K.. "Three dimensional nonlinear dynamic analysis method of underwater line 
structure and its validation", Proceedings of the 10th International Conference on Offshore Mechanics 
and Arctic Engineering (OMAE). ASME, 1991. pp. 87-94. 
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residuals at iteration step / may be expanded in linearized Taylor expansion (in matrix 








form): 
Fe) [Fel [Fe Fa-Fe) [x 
F res = ee + SS X 
Fi E apes OX 
res Jy | res J 
| M + whe, C de K = Ga = GE on; E abo. 
DAt BAt Cx €x E 
E ES ES 
= E q x 
F` Bat öx 
es J 
E a mo 
(3-19) 


- X, 4) and DEM eee X, y) are 


where óx — (Xi P Xie J Ox = (Xa jel 


y»1 
displacement, velocity, and acceleration corrections to be made for the /~ /th iteration at 
time step t+At. An expression for the displacement correction (ôx) may be obtained by 


substituting the last two equations in (3-19) into the first, giving: 


M ele oa K E | -|Z s + Fl «Sm + Ti E =0 
BAt BAt BACCI OX [UM CX ON B CX COD 

(3-20) 
Solving equation (3-20) for 6x and substituting this term back into equation (3-19) gives 
expressions for improved displacement, velocity and acceleration (i.e. displacement, 
velocity and acceleration for iteration step j+ / based upon iteration step J) 


I | 
AA e XA + Ox 
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l 
epo 3 = 
Ne = OX i. (3-21) 
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: CF "— 
The Jacobian matrices, E may be approximated at each iteration step by 
Gx XX 


numerically differentiating the loading vector. 
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With these modifications, the calculation of dynamic equilibrium using the 


Newmark method for the case of highly nonlinear loading could proceed as follows. 


P 


G2 


Assume initial values of displacement, velocity and acceleration (i.e. 


vectors of nodal values). 


For each time step: 


(a) 
(b) 
(c) 
(d) 
(e) 
(f) 


(g) 


Assume initial values of displacement, velocity, and acceleration 
equal to those of the last time step. 

Calculate the loading forces for the current displacement, velocity, 
and acceleration using the chosen physical model. 

Calculate the Jacobian matrices (numerically) for the current 
loading forces. 

Calculate the residuals using equations (3-16) through (3-18). 
Solve for the displacement correction (Ox) using equation (3-20). 
Calculate improved values of displacement, velocity, and 
acceleration using equations (3-21). 

Reiterate through steps (b) through (g) until the residuals are within 


a specified tolerance (i.e. a vector norm tolerance). 


Move onto the next time step, assuming the final displacements, velocities, 


and accelerations of the previous time step. 
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3.3 Definition of dynamic forces. 


3.3.1 Hull structure forces. 

As discussed in section 3.1, the forces associated with the dynamics of the hull 
structure may be represented by an independent formulation from that used to represent 
the internal and external forces. This formulation may involve specific assumptions of 


linearity which may be applied only to the dynamics of the hull structure. 


3.3.1.1 Finite element beam model. 

In representing the whipping response of most conventional ships to low frequency 
bubble pulse loading, it has been found that a simple beam model of the hull structure has 
provided good results’*. This idealization is supported by recognizing that most 
conventional ships (particularly submarines) have length/beam ratios greater than eight. 
Thus, the low frequency vibrational modes occur as beam-like flexural motions, which are 
represented well by simple beam theory. It should be noted, however, that nonlinear hull 
response (i.e. where the assumptions of linear elastic stress-strain relations and small 
rotations cannot be made), makes the utility of the simple beam theory limited. For 
nonlinear hull response, a fully nonlinear, three-dimensional finite element technique would 
be more appropriate. Therefore, for this analysis procedure, the assumptions of linear 
elastic stress-strain relations and small rotations will be made in order to justify the use of 
the simple beam model. 

Because sectional properties vary gradually along the length of the hull, a fre 
element beam model could be used to represent the structural properties of the submarine 
(figure 3-1). The model consists of a specified number of beam elements (equal to the 
number of hull stations specified on the ship's drawings). Most ships have 19 stations 
specified, and thus would be modeled with 19 beam elements. A three-dimensional beam 
model will be used because of the directionality of the loading of the underwater 
explosion, as well as the lack of complete axial symmetry of submarines (i.e. properties 
differ in each plane of symmetry). 

Most methods for structural finite element analysis develop matrices to represent 
the stiffness and mass properties of each finite e/ement, and then combine element matrices 
into global or structure matrices. Modeling of beams using a finite element technique, and 


assignment of the properties of each finite element is discussed in numerous sources 


Hicks. op. cit.. p. 395. 


EM 





(Smith?, Bathe’*, Doyle’, Hughes”, for example). An appropriate finite element analysis 


procedure which could be used to approximate the hull forces associated with submarine 


hull whipping could be outlined as follows: 


D 


D 


dum o. 


v», 


Define the model (coordinate system, nodes, element types, etc.). 
Determine each elements stiffness and mass matrices (k^,m “) and 
transform them to Aul! (or global) coordinates. 

Assemble the Au// or global stiffness and mass matrices (K, M ) by 
superposition of the element matrices. 

Approximate a reasonable Au// damping matrix (C) by applying appropriate 
Rayleigh coefficients, as discussed in sections 2.2.3 and 3.3.1.2 

(C=aM+ BK). 

Combine the hull mass, damping, and stiffness matrices to obtain the total 


hull force as described in section 3.1 (F,,, = MX + Cx + Kx). 


This procedure will be incorporated into an overall computational algorithm and used for 


analyses in the next chapter. 


Elastic beam Lumped 
elements masses/nodes 


Figure 3-1. Structural model for submarine hull whipping. 


Figure 3-2 illustrates a two-dimensional, two-node beam element with vertical 


translation and rotation as the nodal displacements. It should be noted that because a 
submarine should be classified as a relatively thick beam (i.e. the cross-sectional 
dimensions are o/ small compared to the length), the effects of rotary inertia and shear 


deformation at each element should be accounted for in the finite element model. 





73Smith, op. cit. 
“Bathe. op. cit. 
™Dovle. op. cit. 


Hughes. O.F.. Ship Structural Design. SNAME. Jersey City. NJ. 1988. 
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To be rigorous, a three-dimensional beam element could be used which also considers 
additional nodal displacements in the horizontal plane as well as torsion and axial 
displacements (a total of 6 degrees of freedom at each node). Figure 3-3 illustrates a 
general three-dimensional, two-node beam finite element with 6 degrees of freedom at 
each node. It should be noted that axial displacements can generally be ignored for 
submarine whipping due to the lack of axial loading provided by the surrounding fluid. 
Additionally, torsional displacements may also be ignored, although it could be argued that 
torsional loading could become important under some conditions in accounting for 
potential asymmetry associated with the underwater explosion bubble loading on external 
hull appendages such as the sail and planes. Therefore, axial and torsional degrees of 
freedom will be omitted for this analysis, although the entire 3-dimensional beam element 
will be used for illustration in the following discussion. 





Figure 3-3. General three-dimensional, two-node beam finite element. 


Do 





In order to simplify the assembly of the g/oba/ equations of motion for the hull, the 


equations of motion for each beam finite element may be written in the partitioned form: 


ave es C; Ci [JX Ki, Ky tx, _ Jf, 5 
|. Molefe Cs dx: lH Ka Kajlx:] E > 


where subscripts | and 2 refer to nodes | and 2 of the element, and each partition is a 6x6 
submatrix or 6 element subvector. For the general three-dimensional beam element, the 


subvectors for nodal displacement and nodal force would be written in the form: 


a f 
xen f, 
X. f 
la © (3-23) 
O, M, 
0.) Im, | 


where subscript ï refers to the node (1 or 2), x is nodal translation, 9 is nodal rotation, fis 
the (external) force acting on the node, and m is the (external) moment acting on the node. 
It will be shown subsequently that formulation of submatrices for mass and stiffness can be 
readily accomplished. However, the definition of element damping matrices as in equation 
(5-22) 1s generally not possible, thus a formulation of g/oba/ hull damping will be given in 
the next section which will approximate the overall energy dissipation within the hull 
structure. 

The stiffness matrix. The derivation of stiffness properties for a finite element is 
a lengthy process and will not be discussed in detail here. Barltrop and Adams” and 
Hughes”? provide formulation for the element stiffness matrix of a general three- 
dimensional, two-noded beam finite element. The element stiffness matrix can be written 
in a node-oriented partitioned form consistent with equations (3-22) and (3-23). Each 
partitioned submatrix (6x6) has terms i,j which represent the element elastic force in the 
direction of freedom į when a unit displacement is applied to freedom j with all other 
freedoms held at zero displacement. The stiffness submatrices consistent with equation 
(3-22) can be defined as shown ın figure 3-4. 

The mass matrix. The simplest method of assigning mass properties of the 


ship/beam model is to assume that the mass of each section is lumped or concentrated at 


77Barltrop. N.D.P.. and Adams. A.J.. Dynamics of Fixed Marine Structures, Butterworth-Heinemann. 
UK. 1991, Appendix D. 
"8Hughes, op. cit.. pp. 212-217. 
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the nodal points (figure 3-1). This is called the /umped mass formulation. [n this way, the 
lumped masses fill the diagonal terms of the element mass matrix (corresponding to the 
individual nodal degrees of freedom), while the off-diagonal terms are a// zero. The 
matrix (or submatrix) elements would therefore equal the mass of the node for elements 
corresponding to translational degrees of freedom, or the mass moment of inertia of the 
node for those elements corresponding to rotational degiees of freedom. The translation 
of a lumped mass accounts for translational inertia effects, while the rotation of the 
lumped mass accounts for rotary inertia effects. 

The lumped mass formulation has some very appealing properties. Firstly, it is 
easy to associate a physical model with the matrix (each degree of freedom at each node 
having its own mass). Secondly, the entire structure mass matrix is diagonal, leading to 
significantly fewer computations and storage requirements. The finite element mass 
submatrices as defined in equation (3-22), for a lumped mass beam finite element, could be 


written as shown in figure 3-5. 


6l 
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A 1s the cross-sectional area of the element, E is Young's Modulus, L ıs the length of the 


element. I, .I. are the 2nd moments of area of the element cross-section about the z and 


? Zu EN 


y axes, Gis the shear modulus, J is the polar moment of inertia of the element cross- 
section (about the x axis), and A, ,A., are the portions of the cross-sectional area over 


which the y and z-directed shear forces are assumed to act. 


Figure 3-4. Partitioned submatrices for stiffness of a general 3-D beam finite element. 





[0] (no off-diagonal terms for the lumped mass formulation) 


m, is the mass associated with node i 


I.I, are the mass moments of inertia of node i about axes x,v.z. 


Sou Zi 





Figure 3-5. Partitioned submatrices for mass of a general 3-D beam finite element. 


It is also possible to take advantage of the finite element formulation to develop 
what is called a consistent mass formulation for a beam element. In contrast to the 
lumped mass matrix, the components of the consistent mass matrix cannot be associated 
with a physical model. More importantly, after assemblage. the global mass matrix will 
not be diagonal, but will include both diagonal and off-diagonal terms. For this reason. 
the consistent mass matrix formulation requires considerably more computational effort 
than the lumped mass formulation. While the consistent mass matrix formulation is 
generally superior in terms of minimizing errors in calculation of natural frequencies”’, the 
errors associated with the lumped mass formulation significantly decrease as the number of 
nodes used to model the beam increases. Additionally, when external loads are being 
applied directly to the nodes (as will be discussed subsequently), the equations of motion 
behave better when the inertial forces are likewise lumped at the nodes. For these reasons, 


the lumped mass matrix formulation will be used in this submarine whipping model. 


3.3.1.2 Hull damping model. 

Rayleigh damping. As mentioned previously, an approximation for the global 
damping matrix (C) associated with the hull structure can be made by applying 
appropriate Rayleigh coefficients to the global mass and stiffness matrices derived by 


superposition of the element mass and stiffness matrices: 


PDovle. op. cit.. p. 272. 





C=aM+BK (3-24) 
where M and K are the global mass and stiffness matrices, and coefficients « and ß are 
the two damping proportionality constants. As discussed in section 2.2.3, coefficients a 
and p can be calculated if two modal damping ratios (€, and €, ) associated with two 


specific modal frequencies (o , and 0,) are known 


fala 0,0, | ¡QA un D 
D Eod or, buo o. 


The downside to the Ravleigh (or proportional) damping matrix formulation is that 
the contribution of the higher frequency modes of hull flexural response (modal 
frequencies much greater than @ ,) are effectively eliminated by the high damping imposed 
by the selection of the particular coefficients (see figure 2-4). In the case of elastic ship 
whipping, where only the lower few modes of flexural response are generally required to 
adequately describe the motion of the structure, this limitation can be considered 
acceptable as long as the higher frequency used for evaluation of the Rayleigh coefficients 
(0 ,) 1s set among the higher modes expected to contribute significantly to the whipping 
response. Thus, the problem becomes one of calculating the modal (natural) frequencies 
for the lowest modes of hull whipping, applying experimentally determined (or reasonably 
approximated) modal damping ratios for two of the modes, and calculating the appropriate 
Rayleigh coefficients using equation (3-25). As discussed in section 2.2.3, the modal 
(natural) frequencies of a submarine structure can be calculated by solving the generalized 
eigenproblem for the undamped free vibration of the dynamic system: 

Kb, =0,Mb, (3-26) 
where o , is the natural frequency of vibration of the th mode and @, is the mode shape 
vector of the sth mode. Data for experimentally determined (hull) modal damping ratios 
for the fundamental modes of submarine vibrations is lacking in the literature. However, 
some work has been done in the marine industry to estimate hysteretic and structural 
modal damping ratios for steel pile platforms for fundamental-mode flexural vibrations in 
random seas, implying that hysteretic damping in thin steel cylindrical structures is 
typically less than 1% of critical (i.e. € < 0.01) for the fundamental flexural mode. There 
is also some older data available in the literature for modal damping ratios of steel circular 


cylinder ship masts based upon modal "bump" testing?!, suggesting that general steel 


80Cook. M.F.. and Vandiver. J.K.. "Measured and Predicted Dvnamic Response of a Single Pile Platform 
to Random Wave Excitation". Proceedings of the 14th Offshore Technology Conference, OTC 4285, 
1982. pp. 637-645. 

81Smith. op. cit. 
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cylinders have (material and structural) modal damping ratios on the order of 2-3 % of 
critical for the lower flexural modes, although this data seems less reliable because of 
other unknown damping effects on the mast not considered (air and structural connection 
to the ship hull). Although submannes are not precisely the same as thin steel pile 
structures or ship masts, the general hysteretic damping mechanisms are basically the same 
and thus it would seem reasonable to extend this to data for a similar cylindrical steel 
structures to a steel submarine. Thus, it would be reasonable to assume (as a first 
approximation) modal damping ratios (for the hu//) on the order of 1% for the 
fundamental flexural modes of a steel submarine 

Thus, once the modal (natural) frequencies of the hull structure have been 
calculated by equations (5-26), reasonable Rayleigh coefficients could be calculated by 
equations (3-25) assuming modal damping ratios on the order of 0.01. Finally, the hull 
damping matrix to be applied in the finite element calculations can be calculated using the 
global mass and stiffness matrices and equation (5-24). 

This process will be incorporated into the computational algorithm and used in 
analyses in the next chapter. 

An alternative approach. "Dry mode" modal superposition. An alternative to 
developing an explicit damping matrix (C) using Rayleigh damping and solving the 
complete equations of motion, is to use the approximated hull modal damping ratios 
directly and a modal superposition approach. As discussed in section 3.1, each vector 
term of the equations of motion may be represented by an independent formulation (which 
may involve specific assumptions of linearity for that term), so long as all the terms remain 
coupled via a single numerical solution algorithm. Although forces associated with 
external (fluid) loading may be nonlinear, and the equations of motion must therefore be 
solved in the time domain, ıt is still possible to solve for the forces associated with the hull 
structure using a modal superposition technique. This type of approach has been referred 
to as "dry mode" modal superposition??. In a "dry mode" modal superposition approach, 
all forces not associated with the hull are simply termed "loading" forces and thus carried 
on the nght hand side of the equations of motion. 

[MX & CX * Kxj,, = Froosng (DER 
Thus, a coordinate transformation (of the generalized coordinates of the ship huli) may be 
conducted, so long as the transformation of the loading vector (1.€. Fi and Ec, as 


discussed previously) is conducted at each time interval. 


82Bishop. R.E.D.. and Price. W.G.. Hydroelasticitv of Ships. Cambridge University Press. Cambridge. 
UE 1979. 
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This formulation is consistent with the formulation of the Newmark time 
integration method presented in section 3.2 (i.e. equation (3-10)). Because the forces 
associated with the hull are to be modeled using assumptions of linearity (as discussed in 
the previous section), these forces may be redefined using the principles of modal 
superpostion outlined in section 2.2.3. The "dry" mode shapes are first found by solving 
the generalized eigenproblem for the undamped free vibration of the hull structure using 
equation (3-26). The modal equations of motion (equations (2-12)) are then written from 
equations (3-27) by applying the modal transformation, equations (2-10) and (2-13). The 
transformation of the modal force (1.e. Q, in equations (2-13)) must be performed at each 
iteration at each time step of the direct time integration procedure, since the loading 


forces (F ) are generally nonlinear and coupled to the generalized displacements, 


loading 
velocities, and accelerations (i.e. are functions of x,x,X). 

While this approach has some appealing qualities in terms of representation of the 
response of a whipping ship in terms of its fundamental flexural mode shapes, it affords 
little (if any) computational advantage over rigorous computation using the generalized 
coordinates in that two transformations must still be made ar each iteration step of each 
time step. Thus, the fundamental benefit of modal decomposition is lost due to the 
requirement of maintaining the complex (1.e. nonlinear) external loading forces in terms of 
the generalized coordinate system. For this reason, this approach will not be pursued in 


this analysis and all computation will be carried out in terms of generalized coordinates. 


3.3.2 Internal forces. 

In order to analyze the potential contribution of discrete internal structures and 
sloshing liquids on the whipping response of a ship, it is necessary to assume some 
"reasonable" properties of the internal structure and liquid storage tanks as might be 
found on a "real" ship. Appropriate models to account for the dynamics of the internal 


structure and sloshing liquid could be made using simplified geometries and the general 


models discussed in chapter 2. 


3.3.2.1 Model for discrete masses. 
The most straight-forward method of quantitatively determining potential dynamic 


effects of internal structures on the whipping response of a ship 1s to model the internal 
structure with only a few degrees of freedom (i.e. modes of response), and select 
properties of the structure (i.e. mass, stiffness, damping) which are realistic and will make 


the internal structure respond near resonance (response modes no! near resonance would 
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contribute little to the overall response). This amounts to finding a maximum effect 
provided by an internal structure if it were "tuned" to the resonant whipping frequency (or 
the exciting bubble pulse frequency). By using this approach, a quantitative measure of 
the damping effects of internal structures could be obtained for various combinations of 
mass sizes and locations. While mass sizes should represent reasonable ship internal 
structures, the locations of the masses may influence any of the low frequency flexural 
modes of the Aull, and thus a range of potential effects can be predicted. Figure 3-6 
illustrates a notional internal structure connected to the hull through some combination of 


foundation and mounts with known mechanical properties. 





Figure 3-6. Notional internal structure connected to the hull. 


The dynamic vibration absorber. For this analysis, a dynamic vibration absorber 
will be used to model a discrete internal structure. The absorber is considered as a single 
discrete mass, as discussed in section 2.3.1, and connected to the ship hull at one of the 
"nodes" defined in section 3.3.1. Ideally, it would be desirable to permit the absorber to 
respond in all six degrees of freedom (three translational and three rotational). However, 
in order to reasonably limit the computational and modeling effort, the absorber will be 
limited to respond only in the vertical or horizontal directions, as shown in figure 3-7. 
This limitation simplifies the computational algorithm, without unduly restricting the value 
of the model. The absorber model can be thought of as being representative of a resonant 
mode of any discrete internal structure. 

Using the approach developed in section 2.3.1 for a two degree of freedom 
dynamic absorber, the equations of motion for the ship hull, as modified for the forces 


exerted on the hull by the absorber may be written in matrix form: 


| MX t CX = KX] = De ud Due (3-28) 


hull 
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where terms on the left hand sıde are defined by the finite element model for the ship hull 
as described ın section 3.3.1. Displacements and forces are vectors positioned at the hull 
nodes, and hull mass, damping, and stiffness are matrices, constructed for the finite 

element formulatıon. The vector of forces exerted on the hull structure (i.e. at the nodes) 


by the absorber (F, ) is written (ignoring damping associated with the absorber): 


Be ioc m) (3-29) 


where X (capital) refers to the Au// displacement at the node where the absorber acts. and 


bsorber 


x (small case) is the absorber displacement in the direction corresponding to X, as shown 
in figure 3-7. Thus, for a finite element representation of a ship hull with N degrees of 
freedom (i.e. N equations), one additional degree of freedom (i.e. one equation) is 


required for each absorber used. 





Figure 3-7. Dynamic vibration absorber free to respond in 
vertical or horizontal directions. 


Suppression of dynamic response of the whipping hull. As was discussed in 
section 2.3.1, the amplitude of motion of a structure could be decreased if the properties 
of a dynamic vibration absorber were "tuned" to match the resonant response frequency of 
the main structure. If it is desired to consider decrease in amplitude of one of the flexural 
modes of hull whipping by application of a dynamic vibration absorber (the lowest flexural 
mode for example), then the Aul! portion of the two degree of freedom system depicted in 
figure 3-7 could be represented by the modal mass, damping and stiffness for thar 


whipping mode (defined by equations (2-13)). 
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In order to quantify the potential dynamic effects of internal structures on the low 
frequency flexural (whipping) modes of a submarine hull, the following process could be 
implemented to develop the equations of motion for the forces exerted by the absorber 
mass (internal structure) on the hull: 

I: Calculate the mode shapes and natural flexural frequencies of hull whipping 

(€3. ) by solving the finite element eigenproblem as discussed in section 
3.3.1. These represent the dry mode shapes and natural frequencies of 
hull whipping. For a submerged submarine, the mass matrix (M) should be 
modified by adding approximate fluid added mass (M ,) prior to calculating 
the (modal) natural frequencies since the added mass could have a 
significant effect on the natural frequencies. 

2 Calculate the modal mass, damping. and stiffness for the hull whipping 

mode of interest as discussed in section 2.2 2 (equations (2-13)). 


Select a "reasonable" absorber mass or mass of an internal structure (m) 


oJ 


(this could be done for a range of masses). 

4. Calculate the "optimum" absorber stiffness (k), which will "tune" the 
absorber to the resonant frequency of the hull (or the exciting bubble pulse 
frequency): k= m- Q; 

S. Solve for the dynamic response of the structure with the internal absorber 
under specified loading or initial conditions using equations (3-28) and 
(3-29). 

This process will be incorporated into the computational algorithm and used in 


analyses in the next chapter. 


3.3.2.2 Model of liquid sloshing in a rectangular tank. 

Similarly to the case of internal structure, a quantitative measure of the potential 
damping effects provided by internal sloshing liquids may be made by considering the 
forces exerted on the whipping hull by the modal equivalent mechanical system for the 
complex internal sloshing liquid. In effect. only those sloshing modes likely to be excited 
by the whipping of the hull need be considered in the mechanical equivalent system. The 
sloshing modes not excited may be considered as rigid masses, and thus may be considered 
as part of the main hull mass. 

Using the relations developed in section 2.3.2, the equations of motion for the ship 
hull, as modified for the forces exerted on the hull by a sloshing liquid may be written in 


matrix form: 
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[MX +CX+KX] S (3-30) 


hull ^ ^ external 
where terms on the left hand side are defined bv the finite element model for the hull as 
described in section 3.5.1. Displacements and forces are vectors positioned at the hull 
nodes, and hull mass, damping, and stiffness are matrices, constructed for the finite 


element formulation. The vector of forces exerted on the hull structure (i.e at the nodes) 


by the sloshing liquid (F,,,.,, ) can be written (ignoring damping in the tank): 


D --y m, = —mX +) [k,(x, — X)] (3-31) 


where small case letters refer to the mechanical! equivalent dynamic parameters for the 
liquid, and subscript 77 refers to the mode of sloshing. Additionally, the portion of the 
liquid which acts as a rigid body (m,) may be found from equation (2-46) 


m, =M,~ > m, (555/21) 
n=! 


where M, is the total mass of liquid in the tank. 

These relations may be simplified if it is assumed that only the fundamental 
sloshing mode (n = 1) is excited near its resonance (the higher sloshing modes therefore 
may be considered to act as part of the rigid bodv mode, m,). Equation (3-31) and (3-32) 


may be rewritten and combined: 


l 
Frog = —>_m,X, = -(M,X — m,X + m,X,) = (m, - M,)X +[k,(x, - X)] 
n=0 

(9535) 
Thus, for a finite element representation of a ship hull with N degrees of freedom (1.e. N 
equations), one additional degree of freedom (1.e. one equation) is added for each sloshing 
mode of interest and equations (3-30) and (3-35) form a set of N+1 simultaneous 
equations. 

Lateral sloshing in a rectangular tank. Figure 3-8 illustrates a rectangular tank, 
containing a liquid with a free surface, which is oscillating laterally. The potential flow 
solution for the lateral sloshing natural frequencies and hydrodynamic forces are given by 
Vandiver and Mitome* for the case where horizontal displacement of the tank is assumed 
to be harmonic and of the form: 

<a) — Ae. (3-34) 


83 Vandiver and Mitome. op. cit.. pp. 26-28. 
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where A is the amplitude of the harmonic motion of the tank wall. The natural frequencies 


are given for mode n: 





2n-1 
o -aQn- D tanh| : y (3-35) 
a a 


where g 1s the acceleration due to gravity, ais the width of the rectangular tank, and h is 


the mean depth of liquid in the tank. The modal masses are given for each mode n: 


Stanh (2n - |) 1 

Eu ———— 

T n 
a 


As discussed in section 2.3.2, the modal stiffnesses can be computed from the modal mass 


and natural frequency by the expression: 


x 





Figure 3-8. Laterally oscillating rectangular tank with standing waves. 


Suppression of dynamic response of the whipping hull. Similarly to the 
application of the dynamic vibration absorber, the motion of the whipping hull could be 
analyzed by considering any onboard tanks as "liquid slosh" dampers. The main 
differences between sloshing liquids and dynamic vibration absorbers discussed in the 


previous section come about because of the numerous sloshing modes which are not 


A 





excited by the whipping (i.e. the portion of the liquid in the tank which responds as a rigid 
body with the hull, m,). Each sloshing mode which responds may be thought of as an 
individual vibration absorber, but the sloshing modes which do not respond might simply 
be added to the mass of the hull structure. 

Thus, a quantitative feel for the effects of sloshing liquids might be obtained by 
considering the effects of internal structures or absorbers, but also considering the portion 
of the liquids which are not excited, but respond as rigid masses. À more thorough 
discussion of this will be presented in the next chapter. 


3.3.3 External (fluid) forces. 


3.3.3.1 Representation of the hull surface. 

As discussed in section 2.3.1, slender bodies are often represented 
hydrodynamically using a "strip theory" approach since flow around each strip of the hull 
tends to be represented well by considering the flow to be two-dimensional (particularly at 
higher frequencies). Slender submarine hulls, whose length-to-diameter ratios typically 
range from eight to twelve, would therefore be represented well using the hydrodynamic 
"strip theory" approach. 

For purposes of this analysis, a submarine hull will be represented 
hydrodynamically using a number of surface elements ("strips") as shown in figure 
3-9. The centers of pressure of the surface elements will be considered to be located 
coincident with the beam finite element nodes defined in section 3.3.1.1. Thus, the net 
forces exerted on each of the hull sections by the surrounding fluid could be considered to 
act at the nodes of the finite element model. The properties of each surface element (1.e. 
principle dimensions) would be considered constant for the length of that surface element. 

Because two-dimensional flow is assumed around each surface element, the forces 
on each surface element (and therefore on each finite element node) will generally be 
resolved into two Cartesian coordinate components corresponding to the horizontal and 
vertical directions (normal to the ship axis). Force components on each surface element 
will be computed and applied to the beam finite element node. In the case where a 
particular surface element has appendages, additional contributions to the honzontal and 


vertical forces will be computed for the given appendage geometry. 


n» 








Surface element ('strip*) 


Figure 3-9. Hydrodynamic representation of the submarine hull by 
surface elements ("strips"). 


3.3.3.2 Morison relative velocity model. 

As discussed in section 2.3 2, when two dimensional, generally oscillating flow can 
be assumed, and when both fluid inertia and drag forces are important (and wave radiation 
may be ignored), then a Morison's equation may be used to predict hydrodynamic forces 
on a body. It will be assumed for this analysis, that the whipping submarine is sufficiently 
below the free surface such that wave radiation may be ignored. Under this assumption, 
the hydrodynamic forces on each of the surface elements may be predicted using the 
Morison relative velocity formulation given by equation (2-62) and corrected by 


multiplying by the length of the element: 


a " 
Fa 7 [PAC, (à- x) + pAx +>pDC, (u- x)u- v) (3-38) 


where F4, is the total hydrodynamic force on the surface element, p is the fluid density, 


A is the cross-sectional area of the surface element (normal to the direction of flow), D is 
the characteristic diameter of the surface element (normal to the direction of flow), L is 
the length of the surface element, u is the fluid velocity at the location of the center of 
pressure of the surface element (i.e. the "free field" fluid velocity caused by the pulsating 
explosion bubble and considered relatively unaffected by the presence of the hull), and x is 
the displacement of the center of pressure of the surface element (1.e. the finite element 
nodal displacement). C_ and C, are the relative motion fluid inertia and drag 
hydrodynamic coefficients (respectively) and are defined by equation (2-63). 
Determination of hydrodynamic coefficients. For the case of the whipping of a 
submerged submarine, particular care must be taken to select hydrodynamic coefficients 


for analysis which properly represent both the geometry arid the conditions of the relative 





fluid flow around the surface element. Generally, it would be expected that separation of 
flow around the hull would occur during submarine whipping motions, particularly around 
those sections with small characteristic dimensions (or containing components with small 
characteristic dimensions). Determination of the extent to which separation occurs and 
the resulting forces generally requires a complex evaluation of the separation and vortex 
shedding mechanisms. In the "discrete-point" vortex metnod3*, for example, forces on an 
isolated sharp edge mav be predicted my modeling the separating shear lavers which feed 
the growing vortices. The hydrodynamic coefficients are calculated iterativelv over a 
number of cycles by a Fourier averaging process. This type of evaluation. while perhaps 
required for completely rigorous solution of the hydrodynamic forces, 1s beyond the scope 
of this analysis and will not be considered here. 

An initial perspective of the flow fields which could be expected around a 
whipping submarine, and thus of the criteria necessary for initial selection of 
hydrodynamic coefficients for analysis, could be attained by considering past experience 
with ship whipping, and considering typical submarine geometries. As pointed out by 
Hicks®>, the fundamental (i.e. lowest) natural whipping period of most ships is typically 
between 0.3 and 1.0 seconds (fairly well matched with tvpical bubble pulse periods). For 
typical charactenstic dimensions (D) of submarine hull sections (ranges might be between 
0.1 m for small control surfaces, to 15 meters for an entire hull section), characteristic 
frequencies (B) would range from 1 x 10* to 7.5 x 10° (higher whipping modes would 
yield even higher D). These characteristic frequencies correspond to very large Reynolds' 
numbers (Rn) over a large range of Keulegan-Carpenter numbers (KC). 

Extensive experimental studies have been conducted in the marine industry in 
order to obtain values for hydrodynamic coefficients C_, C,, C,, and C}, particularly for 
smooth horizontal cylinders, as discussed in section 2.3.2. Even with the extensive 
studies, a great deal of experimental "scatter" exists in the data base, due to such factors 
as free surface effects, surface roughness, and even the random nature and directionality of 
the ocean environment (for field tests). For these reasons, a great deal of care must be 
taken prior to applying coefficient values from the literature to a design or analysis case. 
Additionally, most experimental studies have been carried out at relatively low frequencies 
(and therefore low Reynolds’ numbers), and thus there ıs limited data which could be 


directly applied to submarıne whipping without some extrapolation. 


84Graham. J.M.R.. "The forces on sharp-edged cvlinders in oscillatory flow at low Keulegan-Carpenter 
numbers". Journal of Fluid Mechanics. 97. 1980. pp. 531-546. 
85 Hicks. op. cit., p. 393. 
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Extensive laboratory research conducted by Sarpkaya3 showed that, for smooth 
circular cylinders at high Reynolds’ numbers (and therefore high ß), the value of C, 


approaches 0.65 and C_ approaches 1.8 (relatively constant over a wide range of KC). 


These coefficients should provide good first approximations for inertia and drag 
coefficients of the circular hull sections (surface elements) for high-frequency whipping. 

Generally, there may be several sharp-edged portions of a submarine surface 
(control surfaces and sail) which may induce significant vortex shedding and therefore 
significant viscous drag on particular hull surface elements Some accounting must 
therefore be made for the additional forces exerted on these surface elements Test data 
on sharp-edged cylinders in high Reynolds’ number oscillatory flow is generally lacking in 
the literature. Some data is presented by Bearman, et. al.’ for flat plates in oscillatory 
flow at Keulegan-Carpenter numbers between | and 10 (relatively low) which might be 
extrapolated to predict that C, would fall in the range of 4-6 and C_ would fall in the 
range of around 1.0 for higher Reynolds’ number (higher frequency) flows. While these 
values are probably not experimentally sound, they should at least provide reasonable first 
order approximations of additional viscous forces exerted on certain hull surface elements 
by sharp-edged control surfaces or sail. 

The modeling of hydrodynamic forces on a whipping submarine using equation 
(3-38) will be included as part of the computational algorithm developed in section 3.4. 
One analysis application will be made in chapter 4 where whipping response of the 
submarine is due only to an initial displacements (i.e. free whipping), and therefore no 
"free field" fluid velocity or acceleration terms (1.e. u and u) will be included, and the 
hydrodynamic forces will be due only to the motion of the hull in otherwise still water. 
Additional analyses will include cases where the submarine is excited by a pulsating 
explosion bubble, and "free field" fluid motion will be included in the hydrodynamic 


model. 


3.3.3.3 Model of fluid loading from a pulsating explosion bubble. 

In analysis cases where it is desirable to couple ship whipping with the loading 
imposed by the pulsation of an explosion bubble, it is necessary to define a time-varying 
formulation for the "free field" fluid velocity and acceleration associated with the pulsating 


86Sarpkava. T.. "In-line and transverse forces on cylinders in oscillating flow at high Reynolds’ numbers", 
Proceedings of the Eight Offshore Technology Conference. OTC 2553. 1976. pp. 95-108. 

87Bearman. P.W.. Downie. M.J.. Graham. J.M.R.. and Obasaju. E.D.. "Forces on cylinders in viscous 
oscillatory flow at low Keulegan-Carpenter numbers", Journal of Fluid Mechanics. 154. 1985. pp. 344- 


345. 





bubble (1.e. u and u). Because derivation of equations for bubble dynamics is a lengthy 
process and not a goal of this analysis, formulation available in the literature will be only 
briefly summarized. A particularly useful formulation for spherically symmetric, pulsating 
and migrating explosion bubbles is presented by Wilkerson®®, which is based upon the 
published theories of Herring??, Friedman??, Taylor?!, Cole?*, and ificks??. This 
formulation includes accounting of effects due to the presence of a free surface and drag 
forces on the migration of the bubble. 

The free field equations of motion for a spherically symmetric, pulsating and 
migrating explosion bubble are described adequately using inviscid. irrotational fluid flow 
theory (1.e. potential flow). Figure 3-10 illustrates an appropriate coordinate svstem which 
encompasses the explosion bubble, its "Image bubble” (to account for the presence of the 
free surface), and an arbitrary point in the fluid. The fluid velocity potential function (o) 
resulting from a pulsating and migrating explosion bubble may be written as the 
superposition of potentials resulting from a simple source and its image (representing 
pulsation of the spherical bubble) and potentials resulting from a dipole and its image 
(representing migration of the spherical bubble) 


e e, e, 
Q - — — —cos0, - — += cos6, (3-39) 
frr epe | E 





where e, and e, are the time-dependent strengths of the sources and dipoles 
(respectively), r, , r. are the radial distances from the center of the bubbles to the arbitrary 
point in the fluid, and 0,,0, are the angles from the vertical (as shown in figure 3-10). 

The first term on the right hand side is the potential due to the source representing the 
pulsation of the explosion bubble. The second term ıs the potential due to the dipole 
representing the migration of the explosion bubble. The last two terms represent the 
potentials due to the effects of the "image bubble" which accounts for the effect of the free 


surface (1.e. satisfies the boundary condition at the free surface). 


58Wilkerson. S.A.. "Elastic Whipping Response of Ships to an Underwater Explosion Loading". Master of 
Science in Engineering Thesis. George Washington University. 1985. 

8?Herring. C.. "Theorv of the Pulsations of the Gas Bubble Produced bv an Underwater Explosion". 
Underwater Explosions Research, Vol. II - The Gas Globe. Office of Naval Research. 1950. 

“Friedman. B.. "Theory of the Underwater Explosion Bubbles". Underwater Explosions Research, Vol. I] 
- The Gas Globe. Office of Naval Research. 1950. 

?1Taylor. G.L. "Vertical Motion of a Spherical Bubble and the Pressure Surrounding It". Underwater 
Explosions Research, Vol. II - The Gas Globe. Office of Naval Research. 1950. 

Cole, R.H.. Undenvater Explosions. Princeton University Press. Princeton. NJ. 1948. 

?3Hicks. A.N.. "The Theory of Explosion-Induced Ship Whipping". Vaval Construction Research 
Establishment, NCRE Report R379. 1972. 
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Figure 3-10. Coordinate system for explosion bubble dynamics. 
Explosion and "image" bubble representation. 


The time-dependent source and dipole strengths are those required to maintain 
continuity of the normal fluid particle velocity at the bubble surface (1.e. the normal fluid 
particle velocity must equal the rate of change of the bubble radius, a, corrected for the 
migration velocity). The required source and dipole strengths are: 


e, =a’a 


ETH 
where v is the upward velocity of the (migrating) bubble and d is the depth of the 


e, ==, : 3 (3-40) 


bubble. Using conservation of energy, the differential equations which govern the 
pulsation and migration of the bubble (i.e. solve for the bubble radius and location) may be 
derived”* The differential equations (in non-dimensionalized form) which describe bubble 


?+Wilkerson, op. cit.. pp. 7-24. 
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pulsation and migration are presented as part of figure 3-11, along with other defined 
constants and non-dimensionalized variables used in the solution of the bubble dynamics. 
The first order differential equations which govern the pulsation and migration of 
the bubble may be integrated using a time-stepping integration procedure such as a Runge- 
Kutta method. It is recommended?” that a variable time-stepping procedure be used 
because of the rapid changes in radius during the phases of bubble minima. With the 
calculated bubble radii and positions (with time), the velocity potential may be calculated 
using equations (3-39) and (5-40), and the "free field" fluid velocity for any arbitrary point 
in the fluid mav be calculated bv differentiating the velocity potential (equation (2-47)) 
Wilkerson’ derived expressions for fluid velocities and accelerations (horizontal and 
vertical components) at an arbitrary point in the fluid as functions of position and source 
and dipole strengths. The resulting expressions (given in Cartesian coordinates of figure 


3-10) are given in figure 5-12. 


5 Hicks. op. cit., p. 159. 
96 Wilkerson. op. cit.. pp. 25-29. 
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Differential equations which govern bubble pulsation and migration: 
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Initial Conditions: a =k it 2 


Non-dimensionalized defined variables: a=a/L (non-dimensional bubble radius) 
==y,/L (non-dimensional bubble depth) 






t=t/T (non-dimensional time) 
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h Ww! 3 \ 
where: L=13.6(W/D,) ° for TNT Ji x] for TNT 





k =(0.0552)D¿ for TNT el 23 for C= 2.25 
Deed +33 dis initial depth of charge W is charge weight of TNT 







Control variables: e (0 for non-migrating; 1 for migrating) 
D (0 for no free-surface effect; 1 for free-surface effect) 





Figure 3-11. Non-dimensionalized variables and differential equations which describe 
the pulsation and migration of an underwater explosion bubble. 
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Figure 3-12. Fluid velocities and accelerations (in relative Cartesian coordinates) for 
given source and dipole strengths and upward bubble velocity. 
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3.4 Formulation of the computational algorithm. 

In order to integrate each of the desired force components presented in the 
previous sections, a comprehensive computational algorithm must be developed. The 
required algorithm must be capable of performing the following: 

1. Compute the time histories of the explosion bubble radius and depth by 

integration of the differential equations governing bubble dynamics using a Runge- 

Kutta time integration procedure. 

2. Compute the "free-field" fluid velocities and accelerations due to the 

pulsating and migrating explosion bubble. These velocities and accelerations must 

be calculated at the locations of each of the finite element nodes used to model the 

dynamics of the hull structure at each iteration at each time step, and must be 
resolved into the vertical and horizontal directions defined by the ship coordinate 
system (a change of coordinates from the bubble coordinate system to the ship 
coordinate system is required). 

3. Compute the external and internal loading on the ship hull using the 

appropriate models developed in the previous sections. 

4. [teratively determine the hull structure displacements, velocities, and 

accelerations using the Newmark method (with Newton-Raphson iterations 

conducted at each time step) as described in section 3.2. Also, iteratively 
determine displacements, velocities, and accelerations of internal masses as 
required. 

A computational algorithm satisfying the above requirements may be generated 
using a suitable computer language, separating the algorithm into sections or subroutines 
appropriate for computation using the chosen language. For this application, the 
algorithm will be implemented on a personal computer using subroutines written for the 
MATLAB® computer program. Appendix A presents MATLAB® subroutines written 
for this analysis, each of which addresses specific portions of the overall algorithm. The 
subroutines are presented as follows: 

a. Subroutine IN.M - Serves as an input file for other subroutines. Data is loaded 
from formatted batch-type files. Arrays and constants used in other subroutines are 
generated. 

b. Subroutine HULL .M - Calculates global mass and stiffness matrices for the hull 
structure, solves the generalized eigenproblem to obtain the "dry mode" natural 
frequencies and mode shapes, calculates an appropriate viscous hull damping matrix for 
the hull using a Rayleigh damping approach, and calculates vectors for modal mass, 


damping, and stiffness. 
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c. Subroutine WETMODE.M - Calculates an approximate fluid added mass 
matrix for the hull using a Morison equation approach (1.e. strip theory), adds the fluid 
added mass matrix to the dry hull mass matrix and solves the generalized eigenproblem to 
obtain the "Wet mode” natural frequencies and mode shapes, and calculates vectors for 
"wet mode" modal mass, damping and stiffness. 

d. Subroutine FREEWHIP.M - Calculates nodal point response (displacements, 
velocities, and accelerations) during free vibration with initial displacements given by 
scaled mode shape vectors for specified modes. Response is calculated using the 
Newmark time integration method with iteration conducted at each time step using a 
Newton-Raphson method. This subroutine is designed primarily to compare damping 
associated with hull mechanisms (i.e. material hysteretic and structural) and damping 
associated with external hydrodynamic forces. 

e. Subroutine BUBBLE.M - Computes and saves the time histories of the 
explosion bubble radius and bubble depth by integration of the differential equations 
governing bubble dynamics. Time histories are computed using a Runge-Kutta integration 
method with a variable time step. The differential equations governing bubble radius are 
called from function routines RADFN.M and RFUNC.M. The differential equations 
governing bubble depth are called from function routines DPTHFN.M and DFUNC.M. 

f. Subroutine BUBWHIP.M - Calculates the nodal point response (displacements, 
velocities, and accelerations) during forced vibration due to a pulsating explosion bubble. 
Response is calculated using the Newmark time integration method with iteration 
conducted at each time step using a Newton-Raphson method. Additionally, local axial 
strains are calculated at specified hull locations using simple beam theory (used for 
comparison of the computational algorithm with model test data). The free-field fluid 
velocities and accelerations due to the pulsating explosion bubble are calculated at each 
iteration when function routine FLUID M is called from within BUBWHIP.M. This 
subroutine is designed primarily to provide a comprehensive computation of explosion- 
induced whipping effects. Additionally, it may be used to quantify damping associated 
with internal structures and to provide for comparison of the computational algorithm with 


model test data. 





4. ANALYSIS OF WHIPPING OF A SUBMERGED SUBMARINE 


EST Introduction. 

[n this chapter, the prediction of the whipping response of a submerged submarine 
Is discussed, including an analysis of specific damping mechanisms. The computational 
algorithm developed in the previous chapter is utilized to predict the response of a 
notional 7400 long ton fast attack submarine (SSN). Additionally, the computational 
algorithm is used to predict the response of the U.S. Navy test platform "Red Snapper”, 
subjected to specified explosive charge weight/standoff combinations for comparison with 
model test data. The latter analysis is presented in non-specitic response units due to 
security concerns with the actual test data used for comparison. 

General characteristics of the notional submarine (SSN) are given in Table 4-1. 
From these general characteristics, specific characteristics (beam element data, hull section 
data, etc.) are formulated and provided as input via subroutine IN.M. The "dry" mode 
natural frequencies and mode shapes are calculated using subroutine HULL.M, and the . 
"wet" mode natural frequencies and mode shapes are calculated using subroutine 
WETMODE.M (t.e. the effects of fluid added mass are included). Table 4-2 gives both 
"dry mode" and "wet mode" natural frequencies computed for the first 12 displacement 
modes. It should be noted that this represents the first 6 displacement modes in each 
plane (i.e. horizontal and vertical), of which the first 2 are considered rigid body modes 
(corresponding to zero natural frequencies), and the remaining 4 modes are considered as 
the fundamental flexural modes. The horizontal "wet" mode shapes are plotted (as 
normalized nodal displacements) and shown in figure 4-1 (rigid body mode shapes) and 


figure 4-2 (fundamental flexural mode shapes). 


Table 4-1. General characteristics of the notional SSN. 


Length (8) 
Displacement/weight (long tons) 7440 


Maximum beam/diameter (ft) 35, 


Hull material HY-80 steel 
Number of hull sections/beam elements 20/19 


Sail. bow planes. stern planes, rudder 

















Table 4-2. Rigid body and fundamental flexural "dry" mode and "wet" mode natural 
frequencies for notional SSN. 


(hz, rad/sec) (hz, rad/sec) 
| Rigid body | 0 | oo 
2 Rigid body | Cd 
3 Rigid body | o | o 
4Rigidbody | | | 
8 Flexural (H) 3.074, 19.315 















Wet mode shape of rigid body mode in horizontal (y=0) plane 
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Figure 4-1. Horizontal rigid body "wet" mode shapes for notional SSN. 
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Wet mode shape of flexural mode 1 in horizontal (y=0) plane 


Wet mode shape of flexural mode 2 in horizontal (y=0) plane 





Wet mode shape of flexural mode 3 in horizontal (yzO) plane 


-1 


Wet mode shape of flexural mode 4 in horizontal (y=0) plane 


Figure 4-2. Horizontal fundamental flexural "wet" mode shapes for notional SSN. 
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4.2 Prediction of "dry" whipping response. Hull damping. 

The "dry" whipping response of the SSN is computed for given initial conditions 
using subroutine FREEWHIP.M (with fluid density set to zero in Morison's equation). 
For free "dry" whipping (as well as free "wet" whipping discussed subsequently), initial 
nodal point displacements are specified as scaled horizontal mode shape vectors. and time 
histories of nodal point displacements, velocities, and accelerations are computed based 
upon these initial displacements. Thus, initial loading on the hull structure is due only to 
internal forces due to structural stiffness. 

For this analysis, initial displacements are set to the first and second horizontal 
fundamental flexural mode shapes (scaled by 2% of the total ship length). Appendix B 
presents plots of horizontal displacements, velocities, and accelerations for nodes 7, 10, 
and 20 for initial displacements (a) scaled from the Ist horizontal fundamental flexural 
mode shape (1.e. mode 6 in Table 4-2), and (b) scaled from the 2nd horizontal fundamental 
flexural mode shape (i.e. mode 8 in Table 4-2). 

As discussed in chapter 3, it has been assumed that the inherent hull damping (i.e. 
due to material hysteresis and structural damping) would reasonably follow experimental 
data for steel cylindrical marine structures and ship masts. Thus, it was assumed that 
modal damping ratios on the order of 1% for hull damping would be sufficient and 
Rayleigh coefficients could be calculated on this basis. Thus, reductions in free response 
amplitude with time (damping) seen in graphs presented in Appendix B are merely 
manifestation of this initial assumption, and no further conclusions may be drawn 
concerning the validity of this initial assumption. 

However, an approximate graphical check could be performed to verify that 
damping in the case of "dry" whipping is indeed represented computationally by the 
assumed 1% modal damping. Indeed, measurement of amplitude log decrements (recall 
equation (1-1)) and calculation of modal damping ratios (recall equation (2-40)) for both 
mode 1 horizontal flexure (represented by node 10 displacement) and mode 2 horizontal 
flexure (represented by node 7 displacement), show that graphically determined modal 
damping ratios are approximately 196. 

4.3 Prediction of "wet" whipping response. Hydrodynamic damping. 

The "wet" whipping response of the SSN is computed for the same given initial 
conditions as the "dry" mode analysis, also using subroutine FREEWHIP.M (with 
inclusion of hydrodynamic forces by Morison's equation). Appendix C presents plots of 
horizontal displacements, velocities, and accelerations for nodes 7, 10, and 20 for initial 
displacements (a) scaled from the 1st horizontal fundamental flexural mode shape, and (b) 


scaled from the 2nd horizontal fundamental flexural mode shape. 
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Figure 4-3 shows comparisons between nodal point displacements for "wet" and 
"dry" cases for initial displacements scaled from mode 1 and mode 2 flexural mode shapes. 
Even though modal damping cannot be directly applied to the nonlinear hydrodynamic 
damping forces, a calculated "equivalent" modal damping might still serve as a measure of 
hydrodynamic damping in comparison to the "dry" case. As was done for "dry" whipping, 
a graphical approximation of "equivalent" modal damping ratios can be made for mode 1 
horizontal flexure and mode 2 horizontal flexure (represented by node 10 and node 7 
displacements, respectively). As expected, the inclusion of external hydrodynamic forces 
(represented by Morison's equation) provide for increased modal damping ratios for both 
flexural modes. Both mode 1 and mode 2 flexure are characterized by "equivalent" modal 
damping ratios of approximately 2.2%. 
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Figure 4-3. Comparison between nodal point displacements for "dry" and "wet" hull 
flexure. Top: mode | flexure, represented by node 10 displacement. 
Bottom: mode 2 flexure, represented by node 7 displacement. 
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4.4 Prediction of whipping response due to explosion bubble pulse loading. 

The previous sections presented analysis based upon an assumption of initial 
displacements set to flexural mode shapes (1.e. free whipping), aimed only at providing 
some quantitative measure of hydrodynamic damping. As discussed in chapter 1, the main 
concern with submarine whipping involves the quasi-periodic bubble pulse loading 
resulting from near to intermediate stand-off explosions It is this quasi-periodic bubble 
pulse loading which can reinforce the free vibration and magnify the whipping response, 
particularly when the period of the bubble pulse is in the vicinity of one of the natural 
periods of hull flexure. Indeed, if the overall dissipative forces (damping forces) are low, 
whipping response magnitudes (and therefore resulting hull damage) could be quite 
substantial. 

To see potential dynamic effects associated with matched bubble pulse (i.e. bubble 
pulse frequency matched to ship's flexural frequency). explosive charge characteristics 
(charge weight, depth, etc.) must first be defined. Using the formulation for bubble 
pulsation and migration discussed in chapter 3, and the known "wet" natural frequencies 
of hull flexure (Table 4-2), explosive charge characteristics could be iteratively 
determined. Table 4-3 specifies explosive charge characteristics used in this analysis to 
match bubble pulse frequency with Ist and 2nd horizontal flexural mode natural 
frequencies (1.e. modes 6 and 8 in Table 4-2). A charge depth of 100 feet was selected, as 
deeper depths require unrealistically large charges to match the Ist fundamental flexural 
frequency (1.37 Hz). Ship location relative to the charge was chosen to promote 
excitation of the mode (both frequency and shape) of interest. 

Figure 4-4 shows plots of bubble radius, depth, and migration velocity for the 
explosive charge weight/depth combination selected to excite mode 2 horizontal whipping 
of the SSN. It can be seen from the plot of radius vs. time that the bubble period ıs nor 
constant with time, but increases as the bubble migrates toward the surface. This is an 
important phenomenon, which limits the ability of the pulsating bubble to excite a 
particular flexural natural frequency for a prolonged period of time (thus, the bubble pulse 
could only be considered to be quasi-periodic in nature). 

The whipping response of the SSN is predicted for the cases when excited by 
explosion bubble pulse characterized by charges specified in Table 4-3. Subroutine 
BUBBLE.M predicts explosion bubble radius and depth time histories for each of the 
charge weight/depth combinations selected. Subroutine BUBWHIP.M predicts the nodal 
point response of the SSN, initially at rest and subject to the incident bubble pulse loading. 
Appendix D presents plots of displacements, velocities, and accelerations for nodes 7, 10, 


and 20 for the charges designed to excite horizontal modes 1 and 2 flexural whipping. 
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Table 4-3. Charge characteristics selected to excite mode | and mode 2 horizontal 
flexural whipping for SSN. 


First horizontal Second horizontal 
flexural en flexural m — 






Charge weight (Ib reach pm | | 1300 | | 
Charge depth (ft) 100 
Bubble pulse frequency (Hz) L35 







Depth of ship (ft) 


[Depthofship(ft) | o% | 80 — 
ship (ft) 
centerline (ft) 
centerline (ft) 


Because of the swiftness with which the (large) bubble migrates to the surface 







from its initial depth for the larger/lower frequency charge (mode | whipping), the 
calculated response time period was limited to the first 2.5 seconds. Despite the short 
time of excitation (only 3 flexure/bubble periods), it can be easily observed that the quasi- 
periodic bubble pulse does indeed magnify the whipping response. For the smaller/higher 
frequency charge (mode 2 whipping), the magnification effect can be more easily observed 


over a larger number of whipping cycles. 
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Figure 4-4. Plots of bubble radius, depth, and migration velocity for the explosive charge 
weight/depth combination selected to excite mode 2 horizontal whipping of the SSN. 
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4.5 Damping from discrete internal structures and sloshing liquids. 

Discrete internal structures. As discussed in chapters 2 and 3, whipping energy 
can be transferred to internal structures within the hull, particularly when the internal 
structure is "tuned" to the bubble pulse frequency (and one of the natural frequencies of 
hull flexure). Using the model developed in chapter 3 (which models the internal structure 
as a discrete mass or dynamic vibration absorber), a brief analysis is conducted to gain 
some insight into potential dynamic effects of these internal structures on the whipping 
response of a hull. 

The analysis discussed in the previous section is extended to include "tuned" 
vibration absorbers (internal masses on springs) where each absorber is "tuned" to the 
matched bubble pulse/whipping frequency. The effects of a "tuned" internal absorber are 
predicted for the case of the bubble pulse matched to the 2nd horizontal flexural 
frequency. The internal absorber is located at node 7 (1.e. "attached" to node 7) and free 
to move in the horizontal direction only, as discussed in section 3.3. Thus, the effect of 
the absorber is on the 2nd horizontal whipping mode. 

The effect of a 200 long ton internal mass, "tuned" to the 2nd horizontal whipping 
mode and "attached" to the hull at node 7, is considered. Appendix E presents plots of 
horizontal displacements. velocities, and accelerations for nodes 7, 10, and 20, as well as 
horizontal displacement, velocity, and acceleration of the absorber mass. Figure 4-5 plots 
horizontal displacement of node 7, with and without the 200 ton internal mass. 

Several interesting phenomenon are evident from the plot and should be briefly 
noted. First, the presence of the absorber has a significant overall effect on reducing the 
magnitude of response. While no specific measure of "damping" can be made, there is 
obviously significant energv transferred into the absorber. Second, the ability of the 
absorber to "absorb" whipping energy is very sensitive to the frequency of excitation. As 
discussed in the previous section, the bubble period (and therefore the frequency of 
whipping excitation) is not constant, but rather increases as the bubble migrates toward 
the surface. By "tuning" the internal absorber to the 2nd horizontal flexural frequency, the 
absorber is not then truly tuned to the excitation frequency (recall that matching bubble 
period to hull flexural frequency was accomplished bv iteration of average bubble periods 
- over the 3 seconds of bubble pulsation). Thus, the usual "internal absorber effect" (i.e. 
reducing vibration amplitude of the main structure to zero) would not occur. Figure 4-6 
plots the horizontal displacement of node 7 and the horizontal displacement of the 
absorber mass. Because the oscillation of node 7 and the oscillation of the absorber mass 
occur at s/ightly different frequencies, a sharp (but temporary) reduction in node 7 
amplitude can be seen around 2.5 seconds, when the absorber shifts from a lag angle to a 
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lead angle (i.e. the node and the mass are temporarily 180° out of phase) What can be 
deduced from this is that internal structures do have potential for significant absorption of 


hull whipping energy, even though bubble pulse frequencies are not truly constant. 
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Figure 4-5. Horizontal displacement of node 7, with and without internal absorber. 
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Figure 4-6. Horizontal displacement of node 7 and 200 LT absorber mass. 
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Liquid sloshing. As discussed ın chapters 2 and 3, the fundamental mechanism by 
which energy is "transferred" to internal liquids through the mechanism of sloshing can be 
thought of much in the same way as that by which it is transferred to internal structures. 
In fact, liquid sloshing may be modeled quite accurately using a mechanical equivalent 
system approach, where each sloshing mode is considered to have its own modal mass, 
damping and stiffness, each responding as its own dynamic system (i.e. its own dynamic 
absorber). Using the formulation for natural sloshing frequency, modal mass, and modal 
stiffness presented in chapter 3 for lateral sloshing in a rectangular tank, the discussion of 
damping through liquid sloshing follows quite nicely behind that of the discrete internal 
structure. Recalling equations (5-35) through (3-37), it 1s clear that the natural frequency, 
modal mass, and modal stiffness of each mode of sloshing are dependent only upon the 
mode number (n), the total mass of liquid in the tank ( M,), and the geometry of the tank 
(1.e. tank depth and width) 

Using the standard rectangular tank formulations. it is a trivial procedure to 
generate plots of natural frequency vs. sloshing mode and sloshing modal mass vs 
sloshing mode. When this is done over a wide range of possible tank configurations, it 
quickly becomes clear that liquid sloshing could not become a Significant contributor to 
the dissipation or absorption of hull whipping energy. Figure 4-7 plots both natural 
sloshing frequency (0 ,) and sloshing modal mass (m,) for a 20' x 20' x 20' rectangular 
tank (a total liquid weight of 228.5 long tons). Two interesting things may be noted for 
this very large tank (much larger than any single, unbaffled tank on any submarine). First, 
the only sloshing modal mass of any substance (i.e. which might effect hull whipping) is 
that of the first/fundamental sloshing mode (i.e. 60 long ton). Second, the natural 
frequency of this fundamental sloshing mode ıs only 2.24 rad/sec (0.35 Hz), well below 
any natural flexural frequency. As other tanks are considered, it is found that the 
fundamental sloshing modes (of all reasonable tank sizes) fall well below 1 Hz. While 
higher sloshing modes may have natural sloshing frequencies in the range of hull whipping 
frequencies, the modal masses of these modes are significantly smaller than would 
contribute to any noticeable reduction in whipping response. Appendix F presents plots of 
natural sloshing frequency vs. sloshing mode and sloshing modal mass vs. sloshing mode 
for several representative rectangular tank geometries. 

From this basic analysis of natural sloshing frequencies and sloshing modal masses, 
it can be justifiably deduced that liquid sloshing would have little impact on the whipping 
response of most ships. Because the significant sloshing modal masses have natural 


frequencies well below hull flexural frequencies, it can also be deduced that the major 





effect of liquids with free surface would only be toward the addition of total ship mass (i.e. 
the significant sloshing modes would act only as rigid bodies). 
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Figure 4-7. Natural sloshing frequency and sloshing modal mass vs. mode number for 
lateral sloshing of a (20' x 20' x 20') rectangular tank. 


4.6 Comparison of computational algorithm with the "Red Snapper" model test. 
The computational algorithm for submarine whipping developed in chapter 3, and 
applied to the whipping of a notional SSN in the previous sections, 1s used to predict the 
response of the U.S. Navy test platform "Red Snapper". Thus, some "verification" of the 
computational algorithm is made with actual model test data. Due to the security 
concerns with the actual test data used for comparison, the following is presented in non- 


specific response units (i.e. the y-axes are unlabeled). 
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The "Red Snapper” model is on the order of 1/4 the size of the notional SSN 
discussed previously (geometrically), and thus natural frequencies of hull flexure are much 
higher than the notional SSN. Whipping response was measured in terms of local axial 
strains (at specific beam element locations) and horizontal and vertical velocities (at 
specific node locations). For comparative purposes, subroutine BUBWHIP.M includes a 
scheme for calculating local axial strains at specified finite beam elements This scheme is 
based upon simple beam theory an provides local axial strains at port and starboard 
extreme fibers, as well as crown and keel extreme fibers for the specified elements (see 
subroutine BUBWHIP.M in Appendix A). 

Three specific model tests, which were conducted on the "Red Snapper” model, 
are used for comparison. One test was designed such that the explosion bubble pulse 
would excite mode 1 horizontal whipping. A second test was designed such that the 
bubble pulse would excite mode | whipping in both horizontal and vertical planes. A third 
test was designed such that the bubble pulse would excite mode 2 honzontal whipping. 
The design of the bubble pulse to match the flexural frequencies and mode shapes was 
carried out in a manner similar to that discussed in section 4.4. In all cases, there where 
no internal structures or liquids, and thus all forces were either hydrodynamic or structural 
in nature. 

General comparisons. Figure 4-8 shows starboard-side fiber strain histones for 
elements 6 and 9 for Test 1 (horizontal mode | whipping), and strains calculated using the 
computational algorithm. In general, it appears that calculated strains match well with 
measured through the 2nd bubble period, but then over predict strains. Figure 4-9 shows 
horizontal velocity histones for nodes 10 and 20 for Test 1, and velocities calculated using 
the computational algonthm. Calculated velocities for this test are only slightly under 
predicted. 

Compansons between measured and calculated strains and velocities for Test 2 
(mode 1 horizontal and vertical whipping) and Test 3 (mode 2 honzontal whipping) are 
presented in Appendix G. In general, the computational algorithm slightly under predicts 
and predicts slightly early for both strains and velocities for Test 2. For Test 3, the 
computational algonthm predicts quite well through the first two bubble periods, but then 


grossly over predicts both strains and velocities. 
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Axial strain vs. time, element 6 
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Figure 4-8. Starboard-side fiber strain histories for elements 6 and 9 for Test 1 of 
the "Red Snapper" model test (horizontal mode | whipping), and strains calculated using 
the computational algonthm 
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Horizontal velocity vs. time, node 10 
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Figure 4-9. Horizontal velocity histories for nodes 10 and 20 for Test 1 of 
the "Red Snapper" model test (horizontal mode 1 whipping), and velocities calculated 
using the computational algorithm 
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Discussion of general comparison. There are several possible explanations for the 
failure of the computational algorithm to properly predict strains and velocities for the 
"Red Snapper” tests. 

a. The bubble model used for formulation of equations governing bubble dynamics 
has been known to be inaccurate past the 2nd bubble period””. Although no satisfactory 
purely physical model exists for long-time bubble dynamics (some more restrictive 
empirical relations are available), this particular bubble model has proven to be fairly 
accurate through the 2nd bubble period. The reasons for the inaccuracy past the 2nd 
bubble period lie in the inability of the model to properly account for dissipation of bubble 
energy during bubble pulsation and migration. Thus, the model believes that the bubble 
maintains more of its original energy throughout the longer time history. This translates 
into greater bubble radii and greater velocities being modeled in the computational 
algorithm (resulting in over predicting hydrodynamic loading past the 2nd bubble period). 
Thus, some explanation for overproduction of response past the 2nd bubble period exists. 

b. The use of the hydrodynamic "line-structure" model in the computational 
algorithm (1.e. hydrodynamic loading being applied at the nodes along the ship centerline, 
vice on the outer hull of the ship), leads to inaccuracies in hydrodynamic loading. This is 
particularly important when the explosive charge is close aboard. In this case, the actual 
hydrodynamic loading is more appropriately seen directly on the outer hull, while the 
computational algorithm believes the loading to be applied at the ship's centerline (i.e. 
calculated "free-field" fluid velocities and accelerations). Thus, an error is introduced due 
to the unaccounted-for distance between the centerline and the outer hull of the ship. For 
large ships, subjected to charges with large standoffs, this geometric inaccuracy Is less 
important. But, for smaller ships subjected to near-standoff charges, this inaccuracy could 
lead to significant error. 

c. There is a general inability of Morison's equation to adequately account for the 
transient-type vortex shedding and convection (viscous drag forces), which is undoubtedly 
important in this application. The original Morison's equation was designed to account for 
hydrodynamic loading for the relatively low-frequency wave loading on pile structures in 
regular sea waves (i.e. quasi-steady). There is no adequate data in the literature to 
support extension of the Morison formulation to transient-type analyses. This is 
particularly important in the justification and selection of the coefficients used in Morison's 
equation. As discussed in chapter 3, coefficients for the computational algorithm were 


merely extrapolated from low frequency data. Additionally, this data was itself obtained 


97Misovec. A.. "Comments on SITE meeting of 5/10/94". (Unpublished meeting minutes). 1994, pp. 1-2. 
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by Fourier time-averaging of measured forces (1.e. quasi-steady). For transient-type 
analyses, such as submarine whipping, a more complex approach would probably be 


appropriate to properly account for the transient vortex shedding and convection. 
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S. CONCLUSIONS 


An analytıcal approach has been presented which predicts the elastic whipping 
response of a submerged submarine to the pulsating bubble of a nearby underwater 
explosion. The approach presented concentrates on mechanisms associated with energy 
dissipation (damping), and includes specific mechanisms associated with hull damping 
(material and structural), internal damping (internal structures and sloshing liquids), and 
external damping (hydrodynamic damping forces in conjunction with traditional inertial 
loading forces). The formulation utilizes a numerical time-domain solution technique to 
directly solve the nonlinear equations of motion. 

While the general mechanisms associated with damping are reasonably well 
understood, application of this understanding to complex high-frequency, non-harmonic 
vibration problems (such as ship whipping due to explosion bubble pulse loading) has not 
been met with much success. For this reason, the approach has been designed to utilize a 
robust solution technique (the time-domain solution), and incorporate physically-based 
models where possible. 

Modal damping ratios for material and structural damping were assumed to be on 
the order of 1 % of critical for fundamental flexural modes, based upon data available 
from steel marine structures. Thus, results for “dry” whipping only reflect this initial 
assumption. With this initial assumption for material and structural damping, “equivalent” 
modal damping ratios for mode 1 and mode 2 whipping (the first two fundamental flexural 
modes) for the case where hydrodynamic damping effects are included (via Morison's 
equation), are computed to be on the order of 2.2 %. 

The effects of internal structures and sloshing liquids on whipping response is 
approximated using a 200 long ton “tuned” mass damper, located to effect mode 2 
whipping. The effects, illustrated in figure 4-5, are significant only when bubble pulse 
frequencies are precisely equal to the “tuned” frequency of the internal structure. The 
effects of internal sloshing liquids are shown to be negligible because of the low natural 
fundamental sloshing frequency associated with larger liquid tanks. 

Several of the selected modeling techniques are shown by comparison with actual 
test data to be questionable. Firstly, the representation of hydrodynamic forces solely by a 
Morison's equation approach fails to fully account for the expected vortex shedding and 
convection forces during the transient vibration. The selection procedure (or more 
appropriately, the calculation procedure) of hydrodynamic coefficients appropriate for ship 
whipping must be better understood and implemented. Mere extrapolation of known data 


based upon harmonic, low frequency motion of small scale cylinders and flat plates is 
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probably inaccurate. A more robust model for hydrodynamic forces, coupling numerical 
simulation with experimental verification should be developed. 

A second modeling technique of questionable success ıs the bubble pulsation and 
migration equations used to produce the free-field fluid velocities and accelerations. 
Because the model does not account for proper dissipation of energy by the pulsating 
explosion bubble, it is believed that the model over predicts fluid loading after the 2nd 
bubble period. A more robust model for fluid loading from the pulsation of the explosion 
bubble, again, coupling numerical simulation with experimental verification should be 
more successful. 

Finally, the general nature of the line-structure model for the ship structural and 
hydrodynamic interaction (i.e. the fluid-structure interaction) introduces geometric errors, 
particularly when the explosion is close aboard. A more robust, 3-dimensional fluid- 
structure interaction solution would be more precise for the solution process. 

Even with several questionable physical models for particular aspects of the 
whipping model, comparison with test data shows that the model provides reasonable 
results, particularly in early-time response. It is believed that the general approach of 
solving the equations of motion directly in the time domain was reasonable and necessary. 
It is in the modeling of each individual component of the forces on the whipping ship 


where future work should be concentrated. 
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% Subroutine IN.M 

% 

% This matlab subroutine ıs written as an INPUT FILE for all other MATLAB 

% subroutines. Data ıs loaded from formatted batch-type input files using 

% the MATLAB "load" command. and arrays and constants are generated for 

% the other MATLAB subroutines. 

% 

% MISCELLANEOUS/DEFAULT VALUES 

% 

% Number of hull stations or finite beam elements 

NE=19: 

% Number of hull sections or finite beam nodes 

NNODES=20: 

% Number of degrees of freedom (+ DOF at each node) 

NDOF=4*NNODES; 

% 

% 

% PROPERTIES OF EACH FINITE BEAM ELEMENT 

load element.dat 

% 

% Element number 

e=element(1.:): 

% Cross-sectional area of element (in2) 

Ae=element(2.:): 

% Length of element (ft) 

Le-element(5.:); 

% Area effective 1n v-directed shear (1n2) 

Asv-element(4.:); 

% Area effective in z-directed shear (in2) 

Asz-element(3,:); 

9o Second moment of area about v-axis (in2ft2) 

Iav-element(6.:). 

% Second moment of area about z-axis (1n2ft2) 

Iaz-element(7.:): 

9$ Y -distance from the neutral axis to the outer fiber of the crown (ft) 

Yc=element(8.:): 

% Y-distance from the neutral anis to the outer fiber of the keel (ft) 
=element(9.:); 

% Z-distance from the centerline to the outer fiber of the port hull (ft) 

Zp=element(10.:): 

% Z-distance from the centerline to the outer fiber of the stbd hull (ft) 

Zs=element(1 1.:): 

% 

% 

% PROPERTIES OF EACH HULL SECTION (NODE) 

load node.dat 

% 

% Node number 

nnode=node( 1.:): 

% Length of each hull section (ft) 

Lh=node(2.:): 

% Weight of hull section (i.e. in air) (Iton) 

wh=node(3.:): 

% mass of hull section (lbf-sec2/ft=slug) 
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mh-7wh*2240./32.174. 

9o Buovancy of hull section (displacement of salt water) (Iton) 
bh=node(+.:): 

% Weight moment of inertia of hull section about y-axis (Iton-ft2) 
Imyh=node(5.:): 

% Mass moment of inertia of hull section about v-axis (slug-ft2) 
Imv-Imvh*2240./52.174; 

“% Weight moment of inertia of hull section about z-axis (Iton-ft2) 
Imzh=node(6.:):; 

“% Mass moment of inertia of hull section about z-axis (slug-ft2) 
Imz-Imzh*2240./32.174. 

% Maximum section beam of main hull perp. to motion in z direction (ft) 
Bz-node(7.:): 

% Maximum section beam of main hull perp. to motion in v direction (ft) 
By=node(8.:): 

% Cross-sectional area of Appended hull (appended portion of each section) 
% perp. to motion in v direction (ft2) 

Aapv-node(9.:); 

% Cross-sectional area of Appended hull (appended portion of each section) 
“% perp. to motion in z direction (ft2) 

Aapz=node(10.:): 

% Width of Appended hull perp. to motion in y direction (ft) 
Dapy=node( 1 1.:): 

% Width of Appended hull perp. to motion in z direction (ft) 
Dapz=node(12.:); 

% Length of Appended hull perp. to motion in y direction (ft) 
Lapv=node(13,:): 

% Length of Appended hull perp. to motion in z direction (ft) 
Lapz-node(14.:): 

% 

% Hydrodynamıc coefficients (non-dımensional) 

% 

% Added mass. inertia. and drag coefficients for each main hull section 
% for motion in vertical (y) and horizontal (z) directions 
Caly=node(135.:); 

Calz=node(16.:): 

Cm ly=node(17.:): 

Cmlz=node(18.:): 

Cdlv=node(19.:): 

Cdlz=node(20.:):; 

% Added mass. inertia. and drag coefficients for appended hull sections 
% (appended portion of each hull section) for motion in y and z directions 
Ca2v=node(21.:): 

Ca2z=node(22.:); 

Cm2v-node(23.:): 

Cm2z=node(2+.:): 

Cd2y=node(25.:). 

Cd2z=node(26.:): 

% 

% 

% MISCELLANEOUS INPUT 

load misc.dat 

% 

% Acceleration due to gravity (ft/sec2) 
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g=misc(1): 

% Density of salt water (Ibf-sec2/ft4=slugs/ft3) 

rho=misc(2): 

% Young's modulus for steel (Ibf/in2) 

E=misc(5): 

% Poisson ratio for steel (non-dimensional) 

Nu=misc(4): 

% Yield stress for steel (Ibf/in2) 

Sigy=misc(5): 

% Shear modulus for steel (lbf/in2) 

G=misc(6): 

% Assumed material (hysteresis) modal viscous damping coefficients for 
% all low frequency (fundamental) flexural "dry" modes (non-dimensional) 
damp=misc(7): 

% Whipping mode to scale as initial condition for use in FREE WHIP.M 
MODE=misc(8). 

% Switch to determine whether or not "wet" whipping is used in FREE WHIP.M 
WET=misc(9): 

y/o Node numbers to plot displacements.velocities.accelerations 
NODE |=misc(10): 

NODE2=misc( 11): 

NODE3=misc(12): 

NODE-4=misc( 15): 

% Element numbers to plot strains 

ELEMENT 1=misc(14); 

ELEMENT 2=misc(15); 

ELEMENT3=misc(16); 

% 

% Input parameters for Newmark time integration algorıthms 
% 

% time step (sec) 

dt=misc(17); 

% gamma parameter (non-dimensional) 

gammap=misc(18); 

% beta parameter (non-dimensional) 

betap=misc(19): 

% time limit (sec) 

tmax=misc(20): 

% vector norm tolerance 

tol=misc(21): 

% 

% Input for internal vibration absorbers 

% 

% "Reasonable" weight of an absorber or internal structure (lton) 
wa=misc(22): 

% mass of absorber or internal structure (slug) 
ma-wa*2240/52.174: 

y/o Node number where absorber is "attached" to the hull 
nabs=misc(23): 

% 

% Input for internal sloshing tank 

% 

y/o "Reasonable" weight of liquid in a (large) internal tank (1ton) 
wl=misc(24): 
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% mass of liquid (slug) 

MI=w1*2240./32.174: 

% "Reasonable" aspect ratio of tank depth/width (h/a) (non-dimensional) 
ha=misc(25); 

% Node number where tank is "attached" to the hull 

ntank=misc(26); 

% 

% Input for explosion bubble dynamics 

% 

% Initial depth of charge (ft) 

Dchg=misc( 27): 

% Charge weight of TNT (Ib) 

Wchg=misc(28): 

% Inıtial locatıon of charge relative to ship hull (in hull coordinate system) 
xchg=misc(29); 

vchg=misc(30): 

zchg=misc(3 1): 

% Initial depth of the ship (ft) 

Dship=misc(32): 


& ccc n dug L 


" dM 
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% SUBROUTINE HULL.M 


% 

% This subroutine performs the following: 

% l. Calculates the global mass and stiffness matrices (K and M) 
% for the 3-D beam finite element formulation (ignoring axial 
% and torsıonal effects). 

% 2. Solves the generalized eigenvalue problem to obtain the 

% "dry mode” eigenvalues (natural freqs squared) and 

% corresponding eigenvectors ("dry mode" shapes). 

% 3. Calculates the Rayleigh coefficients for critical hull 

% (structural and hull hysteresis) damping defined by 

% the Ist and 3rd lowest fundamental flexural modes. 

% +. Calculates the equivalent viscous hull damping matrix for the 
% calculated Ravleigh coefficients. 

% 5. Calculates vectors for modal mass. modal damping. and 

% modal stiffness (for N modes) corresponding to the "drv 
% mode" shape vectors. These are used for optimization 
% of internal absorber and liquid sloshing dampers (tanks) 
% in other subroutines. 

% 6. Provides for validation against the "exact" solution given by 
% Rao (chapter 8) for a continuous uniform free-free beam. 
% The natural frequencies of the finite element solution 

% can be compared to the natural frequencies of the 

% "exact" solution. This may be disabled after validation 
% by placing % before each calculation. 

% 

% (1) 


% Formulate the global (hull) stiffness and mass matrices (K and M). 
% 
% Initialize K and M 
K=zeros(NDOP); 
M=zeros(NDOP); 
% 
% Determine the global stiffness matrix (Ibf/ft for translation or or Ibf-ft for rotation) 
% 
% Determine the element stiffness matrices (KE) 
% 
% Step through each of the 19 elastic beam elements 
for ez I: NE 
% 
Py=(12*E*laz(e))/(G* Asy(e)*Le(e)*2): 
Pz=(12*E*Iay(e))/(G*Asz(e)*Le(e)^2); 
bz=E*laz(e)((1+Py)}*Le(e)^3): 
Dy—E *lav(e)/(( 1+Pz)*Le(e)*3): 
K11=[12*bz 0 0 6*bz*Le(e); 0 12*bv -6*by*Le(e) 0; 
0 -6*by*Le(e) (4+Pz)*by*Le(e)"2 0; 6*bz*Le(e) 0 0 (4+Py)*bz*Le(e)*2]: 
K21=[-12*bz 0 0 -6*bz*Le(e); 0 -12*bv 6*by*Le(e) 0; 
0 -6*by*Le(e) (2-Pz)*bv*Le(e)^2 0; 6*bz*Le(e) 0 0 (2-Py)*bz*Le(e)^2]; 
K12-K2I'*; 
K22=[12*bz 0 0 -6*bz*Le(e). 0 12*by 6*by*Le(e) 0; 
Dosb Tele) ArP7) nv + lele)22 0: -6*bz*Le(e) 0.0 (---Pv)*bz*Le(e)^2]; 


KEZIR TER ID ROT R 22]. 
% 





% Transform element stiffness matrices into their 
% global coordinate matrices 


for 1e=1:8 
for je=1:8 
zierte): 
J=je+4*(e-1):; 
K(1.])) -K(1.)) * KE(1e.je): 
end 
end 
end 
% 


% Once each of the elastic beam elements has been considered and added. 
% the global (hull) stiffness matrix is K. 
% 
% Determine the global mass matrix (Ibf-sec2/ft=slug for translation or slug-ft2 for rotation) 
% 
% Step through each of the 20 nodes 
% 
for n=1:NNODES 
% Step through each of the 4 dof at each node and define the NODE mass matrix 
MN(1,1)=mh(n); 
MN(2.2)=mh(n): 
MN(3,3)=Imy(n): 
MN(4.4)=Imz(n): 
% Transform node mass matrices for each node into their global coordinate matrices 


for in=1:4 
for jn=1:4 
1=1n+4%*(n-1): 
j=jn+4*(n-1): 
Ma.y=May)+MN(«dn,n); 
end 
end 
end 
% 


% Once each of the nodes has been considered. the global 

% mass matrix 1s M 

% 

% (2) 

% Solve the eigenvalue problem to obtain the "dry mode” natural frequencies (square roots of 
% eigenvalues) and corresponding "dry mode” shape vectors (eigenvectors). 

% 

% *** It should be noted that because the hull has 4 rigid body modes, 

% the lowest + eigenvalues will be zero (zero natural frequency). 

% These correspond to the lowest + eigenvectors which are the 

% RIGID BODY mode shapes. Thus, the lowest FLEXURAL mode 

% would be mode n=5. 

% 

% "Phil" is a (NxN) matrix whose columns are "dry mode” shape vectors. 

% "W1" is a diagonal (NxN) matrix whose diagonal elements are (square of the natural frequencies). 
Vo (Note: W is NOT necessarily defined in ascending order). 

% 

[Phil.W]1]-eig(K.M): 

% 

% Define a VECTOR of natural frequencies: 
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for r=1:NDOF 
wl(r)=sgrt(W \(r.r)): 
end 
% 
% Sort natural frequencies and define a vector of natural frequencies given ın ascending order: 
[w2.1]=sort(w1); 
% where "i" is a vector of indices such that w2-w 1(1). 
% 
% Sort the matrix of eigenvectors (mode shapes) as columns in an order corresponding to ascending 
% natural frequencies (defined by index vector "1"): 


for r=1:NDOF 
s=1(r); 
Phi2(:,r)=Phil(:.s): 
end 
% 


% Plot the 2 ngid body mode shapes and the first + flexural mode shapes corresponding to motion in the 
% horizontal plane (for example): 
% Step through the first 12 modes: 


% for m=1:12 

% phi=Phi2(:.m): 

% for n=1:!:NNODES 

% nd(n)=n: 

% ny=1+4*(n-1): 

% nz=2+4*(n-1): 

20 phiy(n)-phi(ny); 

% phiz(n)=phi(nz): 

% end 

% % For rigid body modes (m<=4): 

% if m<=4 

% figure(1) 

% if m<=2, rbm=1: 

% elseif m<=4, rbm>2: 

% end 

% if phiz(1)—0 

% subplot(2,1.rbm), plot(nd.phiz) 
% axis([1 NNODES -1.0 1.0]) 

% title'Drv mode shape of rigid body mode in horizontal (y=0) plane’) 
% end 

% % For flexural modes (m>4): 

% else 

% figure(2) 

% mcount=(m-4); 

% if mcount<=2 

% flexmode=1; 

% elseif mcount<=4 

% flexmode=2: 

% elseif mcount<=6 

% flexmode=3: 

% else. flexmode-4. 

% end 

% if phiz(1)~=0 

% subpiot(4.1.flexmode). plot(nd.phiz) 
% axis([1 NNODES -1.0 1.0]) 

% title('Dry mode shape of flexural mode '.int2str(flexmode).’ in ... 


DS 





% ...horizontal (y=0) plane']) 
% end 

% end 

% end 

% Print "dry" mode shapes 

% figure(1) 


% set(gcf.'PaperPosition',[1.0. 4.5, 6.5. 4.0]) 
% print 

% figure(2) 

% set(gcf.'PaperPosition'.[1.0. 2.0, 6.5, 8.5]) 
% print 

% 

% 

% (3) 


% Calculate the Ravleigh coefficients defined bv the natural frequencies of the Ist and 3rd flexural modes 
?/o (1st and 2nd flexural mode for one of the planes of motion - horizontal or vertical). Ignoring the 4 
% rigid body modes, these would be correspond to modes n=5 and n=7. 

% 

% Define the "dry mode" natural frequencies for modes 5 and 7: 

w5=w2(5): 

w/=w2(7). 

% 

% Find the Rayleigh coefficients using equation (3-16): 

w377[w7 -w5; -1/w7 l1/w5]; 

damping- [damp: damp]: 

pre72*(w5*w7)/(w7^2-w5^2); 

coeff=pre*w57*damping; 

alpha-coeff( 1). 

beta-coeff(2); 

% 

% 

% (4) 

% Calculate the equivalent viscous damping matrix for hull damping using the calculated Rayleigh 
% coefficients: 

% 

C=alpha*M+beta*K: 

% 

% 

% (5) 

% Define the "dry" mode shape vectors and natural frequencies and calculate vectors for modal mass 
% (Mn), modal damping (Cn), and modal stiffness (Kn) corresponding to the "dry mode" mode 

“% shape vectors: 


% 

for n=1:NDOF 
phin=Phi2(:.n); 
phint=phin': 
Phi(:.n)7phin; % "dry" mode shape vectors 
Mn(n)-phint*M*phin; — ?o modal mass vector 
Cn(n)=phint*C*phin; % modal damping vector 
Kn(n)-phint*K*phin; 9 modal stiffness vector 
w(n)-w2(n). % "dry" natural frequencies 

end 

% 

% (6) 
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% Validation of a sample finite element formulation (uniform beam) against the "exact" solution. for a 
% defined axisymmetric (about x-axis) beam arrangement. 

% 

% The "exact" solution is given by Rao (Chapter 8) for a uniform, continuous 

% free-free beam with axisymmetry about the x-axis. This "exact" solution 

% does not include shear effects, so the approximate solution should 

% provide natural frequencies which are slightlv lower than the "exact". 


% 

% Natural frequencies for the "exact" solution for the first 5 transverse flexural modes (modes 1.2.3): 
% 

% Bnl1=4.73: 


% Bnl2=7.85: 

% Bnl3=11.0: 

% rho=mh(1)/Lh(1): 

% Sqroot-sqrt((E*Iav(1))/(rho*(Lh( 1)*NE)^4)): 


% omegal-(Bnl1^2)*Sqroot: 
% omega2-(Bnl2^2)*Sqroot: 
% omega3-(Bnl3^2)*Sqroot: 
% 


% Compare "exact" solution to the (approximate) numerical solution. 

% For approximate solution. rigid body modes are ignored and the Ist. 2nd 

% and 3rd flexural modes for transverse displacement (in each plane (y=0.z=0)) 
% are considered. It should be noted that because of axisymmetry, natural 

“ frequencies and mode shapes in each plane should be the same 

% (1.e. repeated eigenvalues). 


% 

% wl=w(5): 
% w2=w(7): 
% w3=w(9): 
% Compare 

% omegal 

% wl 

% omega2 

% w2 

% omega3 

% w3 


IBI: 





% SUBROUTINE WETMODE.M (called as aMATLAB function) 


% 

% This subroutine performs the following: 

% 1. Calculates the fluid added mass matrix for the submarine 

% using a Morison equation approach (i.e. strip theory). This 

% approach includes added mass effects in the translational 

% degrees of freedom only (rotational fluid added mass effect. 

% are not included). 

% 2. Adds the fluid added mass matrix to the hull mass matrix (subroutine 
% HULL.M) and solves the generalized eigenvalue problem to obtain 
% the "wet mode” eigenvalues (natural frequencies squared) and 
% corresponding "wet mode” eigenvectors ("wet mode" shapes). 
% 3. Calculates vectors for "wet" modal mass. modal damping, and modal 
% stiffness (for N modes) corresponding to the "wet mode" shape 
% vectors. These are used for optimization of internal absorber 

% and liquid sloshing dampers (tanks) in other subroutines. 

% 

% 

% (1) 

% Calculate the fluid added mass matrix (lbf-sec2/ft=slug for translation). 

% 


% Initialize the fluid added mass matrix 
Ma-zeros(NDOF); 
% 
% Step through each of the 20 hull sections (nodes) and define the NODE fluid added mass matrices: 
for n-1:NNODES 
9o Calculate the fluid added mass associated with the main hull (Mal) 
% and the appended hull (Ma2) 
% 
% Step through each of the 4 dof at each node and define the NODE 
“% fluid added mass matrices for the main hull and appended hull 
MNal(1.1)=rho*(pi*By(n)”2/4)*Calv(n)*Lh(n): 
MNa2(1.1)=rho*Aapy(n)*Ca2y(n)*Lapy(n): 
MNal(2.2)-rho*(pi*Bz(n)^2/4)*Calz(n)*Lhí(n): 
MNa2(2.2)=rho* Aapz(n)*Ca2z(n)*Lapz(n): 


MNal(3.3)=0; 
MNa2(3,3)=0: 
MNal(4,4)=0; 
MNa2(4.4)=0: 
% Transform NODE fluid added matrices into their global coordinate matrices 
for in=1:4 
for jn=1:4 
1=1n+4*(n-1): 
j=jn+4*(n-1): 
Ma(i.j)=Ma(i.j)+MNa 1 (in yn)+MNa2(in,jn): 
end 
end 
end 
% 


% Once each of the nodes has been considered. the global fluid added mass matrix is Ma 

% 

% (2) 

% Add the fluid added mass matrix to the hull mass matrix and solve the generalized eigenvalue 
% problem to obtain the "wet mode" eigenvalues (natural frequencies squared) and corresponding 
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% "wet mode" eigenvectors ("wet mode" shapes). 

% 

96 Total mass matrix (Mt): 

Mt=M+Ma: 

% 

% "Phil" ıs a (NxN) matrix whose columns are "wet mode" shape vectors. 

% "W1" 1s a diagonal (NxN) matrix whose diagonal elements are (square of the "wet" natural 
% frequencies). (Note: W is NOT necessarily defined in ascending order). 


% 
|Phil.W1]=e1g(K.Mt): 
% 
% Define a VECTOR of natural frequencies: 
for =1:NDOF 
wl(r)7sqrt(W \(r.r)): 
end 
% 


% Sort natural frequencies and define a vector of "wet" natural frequencies given in ascending order: 
[w2.1]=sort(w 1): 

% where "i" is a vector of indices such that w2=w1(1). 

% 

% Sort the matrix of eigenvectors ("wet" mode shapes) as columns in an order corresponding to ascending 
% "wet" natural frequencies (defined bv index vector "1"): 


for r-1:NDOF 
S-i(r): 
Phi2(:.r)-*Phil (:,s): 
end 
% 


% Plot the 2 rigid body "wet" mode shapes and the first + flexural "wet" mode shapes corresponding to 
% motion in the horizontal plane (for example): 

% Step through the first 12 modes: 

%for m=1:12 


% phi-Phi2(:.m); 

% for n-1:NNODES 

% nd(n)=n: 

% ny=1+4*(n-1); 

% =2+4*(n-]): 

% phiy(n)-phi(ny): 

% phiz(n)-phi(nz); 

% end 

% % For rigid body modes (m<=4): 

% 

% if m<=4 

% figure(1) 

% if m<=2, rbm=1: 

% elseif m<=4, rbm=2: 

% end 

% if any(phiz)—0 

% subplot(2.1.rbm). plot(nd.phiz) 
% axis([1 NNODES -1.0 1.0]) 
% title("Wet mode shape of rigid body mode in horizontal (y=0) plane”) 
% end 

% % For flexural modes (m»4): 

% else 

% figure(2) 
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% mcount=(m-4): 


% if mcount<=2 

% flexmode=1: 

% elseif mcount<=4 

% flexmode=2: 

% elseif mcount<=6 

% flexmode=3: 

% else 

% flexmode-4. 

% end 

% if phiz(1)~=0 

% subplot(+.1.flexmode). plot(nd.phız) 
% axis({1 NNODES -1.0 1.0]) 

% title(["Wet mode shape of flexural mode '.int2str(flexmode).' in 
% ... horizontal (14=0) plane']) 

% end 

% end 

% end 


% Print wet mode shapes 
% figure(1) 


% set(gcf.'PaperPosition'. (1.0. 4.5. 6.5. 4.0]) 
Yo print 

% figure(2) 

% set(gcf,'PaperPosition'.[1.0, 2.0, 6.5. 8.5]) 
% print 

% 

% (3) 


% Define the "wet" mode shape vectors and natural frequencies the first and calculate vectors for "wet 
% modal mass (Mnw), modal damping (Cnw). and modal stiffness (Knw) corresponding to the "wet" 
?o mode shape vectors using equations (2-13): 


% 

for nz 1:NDOF 
phin-Phi2(:.n): 
phint=phin': 
Phiw(:.n)=phin; 9o "wet" mode shape vectors 
Mnw(n)=phint*Mt*phin: % "wet" modal mass vector 
Cnw(n)=phint*C*phin: % "wet" modal damping vector 
Knw(n)=phint*K*phin: % "wet" modal stiffness vector 
ww(n)=w2(n); % "wet" natural frequencies 

end 

% 

% 

% 





% FREEWHIP.M 


% 

% This subroutine performs the following: 

% Calculates the response (nodal point displacements. velocities, 
% and accelerations) of the submarine during a 

% free vibration with initial displacements given bv a 
% scaled mode shape vector for selected mode shapes. 
% Response ıs calculated using the New mark time integration 
% method. with iteration conducted at each time step using a 
% Newton-Raphson method. 

% Utilize the following: 

% - K.C.M.Phi calculated in HULL.M 

% - Hvdrodvnamic coefficients provided in IN.M 

% - Parameters for Newmark algorithm provided in IN.M 

% 

% Define nodes to use for response calculations 

N(1)=NODEI: 

N(2)=NODE2: 

N(3)-NODE3: 

N(43)-NODE4; 

% Define elements to use for axial strain calculations 
Ell=ELEMENTI!; 

EI2=ELEMENT?: 

EIS-ELEMENT3; 

% 


% Set initial and final times for Newmark integration. and limit on number 
% of Newton-Raphson iterations at each time step 
t=0; 
tf-tmax; 
jlimit-50; 
% 
% Initialize the loading forces on the hull 
% Loading forces related to displacement 
Fdisp-zeros(1, NDOF)'; 
% Loading forces related to velocity 
Fvel=zeros(1,NDOF)'; 
% Loading forces related to acceleration 
Facc-zeros( 1. NDOF)'; 
?/o Calculate (constant) static forces (1.e. hvdrostatic/buovancy and weight) 
Fgrav-zeros( Il. NDOF)' 
Fhs-zeros( 1. NDOF)'; 
for n=l: NNODES 
ny=(n-1)*4+1; 
% Hydrostatic/buoyancy force (Ibf) 
Fhs(ny)=bh(n)*2240; 
% Static force of gravity/weight (Ibf) 
Fgrav(ny)=-wh(n)*2240; 
end 
% Total loadıng forces 
Floadt=Fdisp+Fvel+Facc+Fhs+F grav; 
% 
% Define initial conditions (displacements/velocities/accelerations) 
% given by a scaled mode shape vector (for free whipping) 
if MODE==1. xt=(0.02)*sum(Lh)*Phi(:.6): 





else. xt=(0.02)*sum(Lh)*Phi(:.8): 
end 
% 
xdt-zeros( 1. NDOF)'; 
Minv-inv(M); 
xddt=(Minv*(Floadt-C*xdt-K*xt)): 
Yo 
% Perform the Newmark time integration 
Yo 
% Step through time 
Yo 
% Define constants to be used in Newton-Raphson iterations 
cl=(1/(betap*dt”2)): 
c2=(1/(betap*dt)): 
c3=(1/(2*betap)): 
c4=(gammap/(betap*dt)): 
c5=gammap/betap: 
c6=(1-gammap/(2*betap)): 
c7=(1/(Q*betap)-1): 
% 
dx-zeros( 1, NDOF)'; 
% 
for 1=1:300; 
t=t+dt 
% Set initial displacments, velocities,accelerations. and loading forces for 
% the current time step (equal to the final for the last time step) 
X]7xt; 
xdj=xdt: 
xddj=xddt; 
Floadj=Floadt: 
% 
% Perform Newton-Raphson iterations at current time step 
for j=1:jlimit 
% 
% Calculate loading forces for the current iteration 
for n=l:NNODES 
% Loading forces related to velocity 
ny=4*(n-1)+1: 
nz=4*(n-1)+2; 
Fvel(ny)=-0.5*rho*By(n)*Cdlv(n)*xdj(nv)*abs(xdj(nv))*... 
...Lh(n)-0.5*rho* Dapv(n)*Cd2v(n)*xdj(ny)*abs(xdj(ny))*Lapv(n): 
Fvel(nz)=-0.5*rho*Bz(n)*Cd1 z(n)*xdj(nz)*abs(xdj(nz))*... 
...Lh(n)-0.5*rho* Dapz(n)* Cd2z(n)*xdj(nz)*abs(xdj(nz))*Lapz(n): 
% Loading forces related to acceleration 
Facc(ny)--rho*(pi*Bv(n)^2/4)*Calv(n)*Lh(n)*... 
...Xddj(nv)-rho*Aapv(n)*Ca2y(n)*Lapv(n)*xddj(ny): 
Facc(nz)--rho*(pi*Bz(n)^2/4)*Calz(n)*Lh(n)*... 
...Xddj(nz)-rho* Aapz(n)* Ca2z(n)*Lapz(n)*xddj(nz): 
end 
% 
% Total loading force for the current iteration 
Floadj=Fdisp+Fvel+Facc+Fhs+Fgrav: 
% 
% Calculate Jacobian matrices for each of the loading forces 





% (square dıagonal matrices for this formulation) 
dFdxdj 1=zeros(1.NDOF)': 
dFdxddj1=zerost 1.NDOF)': 
dFdxj=zeros(NDOF): 
dFdxdj=zeros(NDOF): 
dFdxddj=zeros( NDOF). 
dxdj=zeros(1,.NDOF)': 
dxddj-zeros( 1. NDOF)': 
for n=1:NNODES 
ny=4*(n-1)+1; 
nz=+*(n-1)+2: 
dxdj(ny)=0.01: 
dxdj(nz)=0.01: 
dxddj(ny)=0.01: 
dxddj(nz)=0.01: 
dFdxdjl(nv)-... 
...(-0.5*rho*Bv(n)*Cdlv(n)*(xdj(ny)*dxdj(ny))*abs(xdj(ny)*dxdj(ny))*... 
...Lh(n)- 
0.5*rho*Dapv(n)*Cd2v(n)*(xdj(nv)*dxdj(ny))*abs(xdj(nv)*dxdj(ny))*... 
... Lapy(n)-Fvel(ny))/dxdj(ny): 
dFdxdj1(nz)=... 
...(-.5* rho*Bz(n)*Cd1z(n)*(xdj(nz)+dxdj(nz))*abs(xdj(nz)+dxdj(nz))*... 
...Lh(n)- 
0.5*rho*Dapz(n)*Cd2z(n)*(xdj(nz)*dxdj(nz))*abs(xdj(nz)*dxdj(nz))*... 
...Lapz(n)-Fvel(nz))/dxdj(nz ); 
dFdxddjl(nv)-(-rho*(pi*By(n)^2/4)*Calv(n)*Lhí(n)*... 
...(xddj(nv)*dxddj(ny))-rho* Aapv(n)*Ca2y(n)*Lapy(n)*... 
..(xddj(ny)+dxddj(ny))-Facc(ny))/dxddj(ny); 
dFdxddj1(nz)=(-rho*(pi*Bz(n)*2/4)*Calz(n)*Lh(n)*... 
...(xddj(nz)*dxddj(nz))-rho* Aapz(n)*Ca2z(n)*Lapz(n)*... 
..(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz): 


end 

for n=1:NDOF 
dFdxdj(n.n)=dFdxdj1(n): 
dFdxddj(n.n)-dFdxddj l (n); 

end 

% 


% Calculate residuals 

f1j=M*(cl *(xj-xt)-c2*xdt-c3 *xddt)+C*(c4*(xj-xt)-C5 *xdt+c6*xddt*dt)+... 
...K*(xj-xt)-Floadj+Floadt: 

f2j-xdj-c4*(xj-xt)* (c5- 1)*xdt-có6*dt*xddt; 

f3j7xddj-cl*(xj-xt)* c2 *xdtc7* xddt: 

% 


% Solve for the displacement correction (dx) 

% 
A1=(c1*M+c4*C+K-<1*dFdxddj-c4*dFdxdj-dFdxj): 
Alinv=inv(A 1): 

A2=-f1j-dFdxdj* f2j-dFdxddj*f3): 

dx=(Alinv*A2); 

% 

% Calculate the next iteration's displacements. velocities, and accelerations 
% 

Xj=Xjtdx, 

xdj=xdj+(c4*dx)-f2): 





xddj=xddjt+(c1 *dx)-f3): 
% 
% Check to see if all Norms of residuals are less than tolerance for 
% the current iteration 
% (If so. go to the next time step) 
Normf1j-norm(f1]): 
Normf2j=norm(f2)): 
Normf3j=norm(f3)): 
if Normflj<=tol & Normf2j<=tol & Normf3j<=tol 
% 
% Save current time, displacements, velocities. 
% and accelerations (at specified nodes) for later plotting 
time(1)=t: 
% Displacements/velocities/accelerations 
for n=1:4 
Dh(n)=(N(n)-1)*4+2: 
end 
X1(i)=xj(Dh(1)): 
Xd1(i)=xdj(Dh(1)): 
Xdd 1(1)=xddj(Dh(1)): 
X2(i)=xj(Dh(2)): 
Xd2(i)=xdj(Dh(2)); 
Xdd2(i)=xddj(Dh(2)): 
X3(i)=xj(Dh(3)): 
Xd3(1)=xdj(Dh(3)): 
Xdd3(1)=xddj(Dh(3)): 
X4(i)=x)(Dh(4)); 
Xd4(i)=xdj(Dh(4)). 
Xdd4(i)=xddj(Dh(4)): 
% 
Dh19=(19-1)*4+2; 
Dh18=(18-1)*4+2; 
Fvel18(1=Fvel(Dh18): 
Facc18(1)=Facc(Dh18); 
Fload18(i)=Floadj(Dh18): 
Fvel19(1)=Fvel(Dh1 9). 
Facc19(1)=Facc(Dh 19). 
Fload 19(1)=Floadj(Dh1 9): 
% 
% Save the current time step displacements, velocities, and 
% accelerations, and forces for use in the next time step 
Xt=Ny: 
xdt=xdj: 
xddt=xddj: 
Floadt=Fload): 
break. end 
end 
if j==jlimit, break, end 
end 
if j==jlimit. break, end 
if t>=tf. break, end 
end 
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% Subroutine BUBBLE.M 


% 

% This subroutine performs the following: 

% 1. Computes the time histories of the explosion bubble radius (rad) 

% and bubble depth (depth) by integration of the differential 

% equations governing bubble dynamics. 

% Computation is carried out for specified bubble initial conditions 
% (i.e. initial depth of the charge (Dchg) and weight of the 

% charge (Wchg)). 

% 

% 


function[TM.BR.BD.BRD.BDD.BRDD.BDDD.DEL.btime.brad.bdpth.bvel.L1.T1]=... 
...bubble(Dchg. Wchg.tmax) 

% 

% Define constants 

% 

% Drag coeff. for bubble 

Cd=2.25; 

% Initial bubble head (ft) 
DO=Dchg+33; 

% Adiabatic gas exponent for TNT 
gl=1.25; 

% Constant for TNT 
k12(0.0552)*D0^0.25; 

% Length scale constant 
L1=13.67*(Wchg/D0)(1/3): 

% Time scale constant 
T122.94*(Wchg^(1/3//D0^(5/6)); 

“% Misc. constants 

cl233/L1; 

c2=D0/L 1; 

c3=(g1-1)*k1; 

c4=3*g1+1: 

% 

% Define initial and final conditions 
% 

% Non-dim. bubble radius 

x0=k 14(4/3)*(1+(4*k1%4)/3); 

% Non-dim. bubble depth 

NU=SC2: 

% Rate of change of bubble radius 
xd020; 

% Rate of change of bubble depth 
vd0-0: 

% Non-dim. start time (i.e. 0 seconds) 
TO=0.0; 

% Non-dim. end time (1.e. tmax seconds) 
Tf=tman/T 1. 

% Non-dimensional time step 
dT-7(5.0*x0^4)/Tf. 

% Integration time step parameter 
tsp=10.: 

% 

% Perform Runge-Kutta time iteration 





% 


T=TO0: 
X-x0; 
SEND, 
xd=xd0: 
vd=yd0: 
% 
M=3000: % defined maximum number of time steps (for control) 
% 
for m=1:M 
% Call the differential eqn. for bubble radius 
[x.xd]=radfn(dT.x.xd.y.yd.c1.c2.c3.c4.Cd): 
xdd=-1.5/(1-x/(2*(y-c1)))* ((xd^2/x)*(1-2*x/(3*(y-c 1)))-vd^2/(6*x)+y/(x*c2)-c3/(x^c4)+... 
..(x/(4*(y-c1)^2))*(xd*yd/3-x/c2+Cd*yd^2/4)): 
% Call the differential eqn. for bubble depth 
[v.yd]=dpthfn(dT.x.xd.v.vd.xdd.c1.c2.c3.c4.Cd): 
ydd--3*(1/c2*xd*vd/x-(Cd*vd^2)/(4* x)*(x/(4*(y-c1)^2)) *(3*xd^2*x* xdd)): 
% 
% Save the current bubble characteristics for plotting and later use 
% 
TM(m)=T: % non-dim. time 
BR(m)=x: % non-dim. bubble radius 
BD(m)=y; % non-dim. bubble depth 
BRD(m)=xd: % non-dim. rate of change of bubble radius 
BDD(m)=yd: % non-dim. rate of change of bubble depth (bubble velocity) 
BRDD(m)=xdd: 
BDDD(m)=vdd. 
DEL(m)=(y-cl): % non-dim. suface pressure 
% 
btime(m)=T*T1; % dimensional time 
brad(m)=x*L1; % dimensional bubble radius 
bdpth(m)=(v-c1)*L1: % dimensional bubble depth 
bvel(m)=vd*L1/T1: % dimensional bubble velocity 
% 
% Go to next time increment 
x3=x 83: 
x3d73.0*(x^2)*xd; 
x3dd=3.0*x*(x*xdd+2.0*xd"2); 
dT=2.0/(abs(x3dd)*3.0*tsp): 
if AT<0.001, dT=0.001; 
elseif dT>0.025, dT=0.025: 
end 
T=T+dT: 
if T>Tf. break 
end 
end 
% 
% Plot the results 
% 
figure(1) 


subplot(3.1.1). plot(btime.brad) 
title('Bubble radius vs. time”) 
xlabel("Time (sec)') 
vlabel('Bubble radius (ft)') 
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subplot(3.1.2). plot(btime.-bdpth) 
tide('Bubble depth vs. time”) 

xlabel("Time (sec)') 

vlabel('Bubble depth (- ft)') 
subplot(3.1.3). plot(btime. -bvel) 
title('Upward bubble velocity vs. time’) 
xlabel("Time (sec)') 

vlabel('Upward bubble velocity (ft/s)') 
set(gcf.'PaperPosition'.[1.0, 2.0, 6.5, 8.5)) 
% print 





% RADFN.M 

% 

% Function M-file for calculatıon of the bubble radius 

% (1.e. contains the bubble radius governing equation for use in BUBBLE.M) 
% 

function [x.xd]=radfn(dT.x.xd.y.yd.c1,c2.c3.c4.Cd) 

ak 1 =dT*rfunc(x.xd.v.vd.c1.c2.c3.c4.Cd): 
ak2-dT*rfunc(x*dT*xd/2.0.xd-*ak1/2.0.v.vd.c1.c2.c3.c4.Cd). 
ak3-dT*rfunc(x*dT*xd/2.0*-dT*ak1/4.0.xd*ak2/2.0.v.vd.c1.c2.c3.c4.Cd): 
ak4=dT *rfunc(x+dT*xd+dT *ak2/2.0.xd+ak3.v.vd.cl.c2.c3.c4.Cd): 
x=x+dT*xd+dT *(ak 1 +ak2+ak3)/6.0: 
xd=xd+(ak1+2.0*ak2+2.0*ak3+ak4)/6.0: 


% RFUNC.M 

% 

function ak=rfunc(x.xd.v,yd.c1.c2.c3.c4.Cd) 

ak=-1.5/(1-x/(2*(y-c1)))*((xd“2/x)*(1-2*x/(3 ¥(v-c1)))-vd*2/(6*x)+W/(x*C2)-€3/(x%C4)+(X/(4 *(y- 
ely 2))* (xd *yd/3-x/c2+Cd*yd%2/4)). 
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% DPTHFN.M 

% 

% Function M-file for calculation of the bubble depth 

% (i.e. contains the bubble depth governing equation for use in BUBBLE. M) 
% 

function [v.vd]=dpthfn(dT.x.xd.v.vd.xdd.cl.c2.c3.c4.Cd) 
akl-dT*dfunc(x.xd.v.vd.xdd.cl.c2.c3.c4.Cd): 
ak2-dT*dfunc(x.xd.v*dT*vd/2.0.vd*ak1/2.0.xdd.c1.c2.c3.c4. Cd): 
ak3-dT*dfunc(x.xd.v*-dT*vd/2.0*dT*ak1/4.0.vd*ak2/2.0.xdd.c1.c2.c3.c4. Cd). 
ak4-dT*dfunc(x.xd.v-*dT*vd-*dT*ak2/2.0.vd*ak3.xdd.c1.c2.c3.c4. Cd). 
v=vt+dT *yd+dT*(ak | +ak2+ak3)/6.0. 
vd=yd+(ak1+2.0*ak2+2.0*ak3+ak4)/6.0; 


% DFUNC.M 

% 

% 

% 

function ak=dfunc(x.xd.y.yd.xdd.cl.c2.c3.c4.Cd) 
ak--3*(1/c2*xd*vd/x-(Cd*vd^2)/(4*x)*(N/(4*(v-c1)^2))*(3*xd^2*x*xdd)): 





% BUBWHIP.M 


% 

% 

% This subroutine performs the following: 

% l. Calculates the response (nodal point displacements. velocities. accelerations) 
% of the submarine during a vibration with initial displacements of zero 
% and subjected to fluid loadıng produced by a pulsatıng explosion bubble. 
% Response is calculated using the Newmark time integration algorithm, 
% with iteration conducted at each time step using a Newton-Raphson method. 
% 2. Calculates the local axial strains at specified hull locations using simple 

% beam theory (used for comparison with model test data). 

% 

% 


% Utilize the following: 

% -K.C.M calculated in HULL.M 

% -Hydrodynamic coefficients provided in IN.M 

% -Parameters for Newmark algorithm provided in IN.M 

% -Bubble characteristics (radius. depth, etc.) calculated in BUBBLE.M 

% 

% Define nodes to use for response calculations 

N(1)=NODEI; 

N(2)=NODE?: 

N(3)-NODE3; 

N(4)=NODE4; 

% Define hull elements to use for axial strain calculations 

EIIZELEMENTI: 

EI2=ELEMENT?: 

EI3=ELEMENT3; 

% 

% Define the horizontal (hull) degree of freedom where internal structure/absorber acts 

nnabs=(nabs-1)*4+2; 

% 

% Define parameters for internal structure/absorber 

ka-ma*(ww(8)^2); 

mabs-zeros( 1. NDOF)'; 

kabs-zeros( 1. NDOF)' 

mabs(nnabs)=ma; 

kabs(nnabs)=ka: 

% 

% Set the initial and final times (sec) for Newmark integration, and limit 

% on the number of Newton-Raphson iterations 

=(): 

tf-tmax: 

jlimit=40; 

% 

% Initialize the loading forces on the hull 

% Loading forces related to displacement 
Fdisp=zeros(1.NDOF”)': 

% Loading forces related to velocity 
Fvel=zeros(1.NDOF)’: 

% Loading forces related to acceleration 
Facc=zeros(1,.NDOF)'. 

% Calculate (constant) static forces (i.e. hydrostatic/ouovancy and weight) 

Fgrav-zeros( 1. NDOF)'; 
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Fhs-zeros( l.NDOF)': 
for n-1:NNODES 
ny=(n-1)*4+1; 
% Hydrostatic/buoyancy force (Ibf) 
Fhs(nv)-bh(n)*2240: 
% Static force of gravitv/weight (Ibf) 
Fgrav(ny)=-wh(n)*2240: 
end 
% Total initial loading forces 
Floadt=Fdisp+Fvel+Facc+Fhs+Fgrav; 
% 
% Define initial conditions (displacements/velocities/accelerations) 
% Hull 
xt=zeros(1,NDOF)': 
xdt=zeros(1.NDOF)': 
xddt-zeros( 1L NDOF)': 
% Internal absorber 
xat-0; 
xadt-0; 
xaddt-0; 
% 
% Perform the Newmark time integration 
% 
% Step through time 
% 
% Define constants to be used in Newton-Raphson iterations 
cl=(1/(betap*dt”2)): 
c2=(1/(betap*dt)): 
c3=(1/(2*betap)): 
c4=(gammap/(betap*dt)): 
c5=gammap/betap: 
c6-(1-gammap/(2*betap)): 
c7=(1/(2*betap)-1): 
% 
dx=zeros(1,.NDOF)': 
% 
for 1=1:199; 
t=t+dt 
% 
% Set initial displacments, velocities,accelerations, and loading forces for 
% the current time step (equal to the final for the last time step) 
xj7xt: 
xdj=xdt; 
xddj=xddt; 
Floadj=Floadt. 
xaj=xat; 
xadj=xadt: 
xadds=xaddt: 
% 
% Perform Newton-Raphson iterations at current time step 
for j- 1 limit 
% 
xaddj=xadds: 
% Calculate current iteration absorber displacement/Velocity 





xaj=xat+xadt*dt+xaddt*(dt”2/2)+betap*(xaddj-xaddt)*dt"2; 
xadj=xadt+xaddt*dt+gammap*(xaddj-xaddt)*dt: 

% 

% Calculate the loading forces for the current iteration 

% 

% Call function "fluid" to provide free field fluid velocities and 

“% accelerations at each node (at the current iteration node location) 


[uv.uz.uyd.uzd]-fluid(xj.xchg.vchg.zchg.Dship.t. TM.BR.BD.BRD.BDD.BRDD.BDDD.... 
...DEL.L I.T I.btime.bdpth. NDOF.NNODES.Lh); 

% 

for n=1:NNODES 
nv=4*(n-1)+1: 
nz=4*(n-1)+2: 

% Loading forces related to displacement 
Fdisp(nz)=kabs(nz)*(xaj-xj(nz)). 

% Loading forces related to velocity 
Fvel(ny)=0.5*rho*By(n)*Cdly(n)*(uv(n)-xdj(ny))*... 
...abs(uy(n)-xdj(ny)) *Lh(n)*0.5*rho*Dapv(n)*Cd2y(n)*(uv(n)- 
...Xdj(ny))*abs(uy(n)-xdj(ny))*Lapy(n): 
Fvel(nz)=0.5*rho*Bz(n)*Cdlz(n)*(uz(n)-xdj(nz))*... 
...abs(uz(n)-xdj(nz))*Lh(n)*0.5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 
...Xdj(nz))*abs(uz(n)-xdj(nz))*Lapz(n); 

% Loading forces related to acceleration 
Facc(ny)-7rho*(pi*Bv(n)^2/4)*Cmly(n)*Lh(n)*(uyd(n)- 
...Xddj(ny))* rho* Aapy(n)*Cm2y(n)*Lapy(n)*(uvd(n)- 

...Xddj(ny))*rho*(pi*Bv(n)^2/4)*Lh(n)*xddj(ny)*rho* Aapy(n)*Lapy(n)*xddj(nv). 
Facc(nz)-rho*(pi*Bz(n)^2/4)*Cmlz(n)*Lh(n)*(uzd(n)- 
...xddj(nz))*rho* Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 

...xddj(nz))* rho*(pi*Bz(n)^2/4)*Lh(n)*xddj(nz)*rho* Aapz(n)*Lapz(n)* xddj(nz). 

end 

% 

% Total loading force for the current iteration 

Floadj=Fdisp+Fvel+Facc+Fhs+F grav; 

% 

% Calculate Jacobian matrices for each of the loading forces 

% (square diagonal matrices for this formulation) 

% 

dF dxj 1=zeros(1.NDOF)’: 

dFdxdjl-zeros(1, NDOF)': 

dFdxddjl-zeros(1. NDOF)': 

dFdxj-zeros(NDOF); 

dFdxdj-zeros(NDOF): 

dFdxddj=zeros(NDOF): 

dxj=zeros(1.NDOF)’: 

dxdj=zeros(1,NDOF)’: 

dxddj=zeros(1.NDOF)': 

for n=1:NNODES 
ny=4*(n-1)+1; 
nz-4*(n-1)*2; 
dxj(nz)=0.01; 
dxdj(ny)=0.01; 
dxdj(nz)=0.01: 
dxddj(nv)-0.01; 





dxddj(nz)=0.01: 

dFdxj l(nz)-(kabs(nz)*(xaj-(xj(nz)*dxj( nz)))-Fdisp( nz))/dxj( nz): 

dFdxdjl(nv)-(0.5*rho*Bv(n)*Cdlv(n)*(uv(n)-(xdj( ny)*dxdj(ny))) *abs(uv(n)- 

...(xdj(ny)*dxdj(nv))) *Lh(n)*0.5*rho*Dapv(n)*Cd2v(n)*(uv(n)- 

...(xdj(ny)*dxdj(nv)))*abs(uv(n)-(xdj(nv)*dxdj( nv)))*Lapy(n)- 

...Fvel(ny)Ydxdj(ny): 

dFdxdjl(nz)-(0.5*rho*Bz(n)*Cdlz(n)*(uz(n)-(xaj(nz)*dxdj(nz))) *abs(uz(n)- 

...(xdj(nz)*dxdj(nz))) *Lh(n)70.5*rho*Dapz(n)*Cd2z(n)*(uz(n)- 

...(xdj(nz)*dxdj(nz)))*abs(uz(n)-(xdj(nz)*dxdi(nz))) *Lapz(n)- 

...Fvel(nz))/dxdj(n2): 

dFdxddjl(nv)-(rho*(pi*Bv(n)^2/4)*Cml v(n)*Lh(n)*(uvd(n)- 

...(xddj(nv)*dxddj(ny)))*rho* Aapy(n)*Cm2v(n)*Lapv(n)*(uvd(n)- 
...(xddj(nv)*dxddj(nv)))*rho*(pi*By(n)^2/4)*Lh(n)*( xddj(nv)*dxddj( nv))*rho* Aapy(n)* .. 

... Lapy(n)*(xddj(nv)+dxddj(nv))-Facc(ny))/dxddj(ny): 

dFdxddjl(nz)-(rho*(pi*Bz(n)^2/4)*Cmlz(n)*Lh(n)*(uzd(n)- 

...(xddj(nz)+dxddj(nz)))+rho* Aapz(n)*Cm2z(n)*Lapz(n)*(uzd(n)- 
...(xddj(nz)*dxddj(nz)))*rho*(pi*Bz(n)^2/4)*Lh(n)*(xddj(nz)*dxddj(nz))*rho* Aapz(n)* ... 

... Lapz(n)*(xddj(nz)+dxddj(nz))-Facc(nz))/dxddj(nz). 


end 

for n=1:NDOF 
dFdxj(n.n)=dFdxj1(n): 
dFdxdj(n.n)-dFdxdj (n); 
dFdxddj(n.n)-dFdxdd]l (n); 

end 

% 


% Calculate residuals 

f1j=M*(c1 *(xj-xt)-c2 *xdt-c3 *xddt)+C *(c4 *(xj-Nt)-c5 *xdt+c6*xddt*dt)+... 
...K*(xj-xt)-Floadj+Floadt: 

f2j=xdj-c4 *(xj-xt)+(c5-1)*xdt-c6 *dt*xddt; 
Bj=xddj-c1*(xj-xt)+c2*xdt+c7*xddt: 

% 

% Solve for the displacement correction (dx) 

% 

A1=(c1*M+c4*C+K<1*dFdxddj-c4*dFdxdj-dFdxj): 

A linv=inv(A 1); 

A2=-flj-dFdxdj*f2j-dFdxddj*f3j: 

dx-(A linv* A2): 

% 

% Calculate the next iteration's displacements. velocities, and accelerations 
% 

x]=xj+dx; 

xdj=xdj+(c4 *dx)-f2): 

xddj=xddj+(cl *dx)-f3]; 

% 

% Required acceleration of absorber 

if any(mabs)~=0, xadds=-Fdisp(nnabs)/mabs(nnabs). 

else, xadds=0; 

end 

% 

% Check to see if all Norms of residuals are less than tolerance for 
% the current iteration 

% (lf so. go to the next time step) 

Normflj=norm(flj): 

Normf2j-norm(f2j). 





Normf3j=norm(ß3)): 
Normabs=norm(xadds-xadd)): 
if Normf1lj<=tol & Normf2j<=tol & Normf3j<=tol & Normabs<=tol 
% 
% Calculate axial stress at specified hull stations 
% (using simple beam theory) 
% 
qy=zeros(1.NE)': 
qz-zeros( 1. NE)'; 
Vy-zeros(1.NE)': 
Vz-zeros( l.NE)': 
Myl1-zeros(1.NE)': 
Mzl-zeros( 1l. NE); 
My=zeros(1.NE)': 
Mz-zeros(1.NE)'; 
for e- 1: NE 
% Net vertical or horizontal force at element e 
ey=(e-1)*4+1: 
ez=(e-1)*4+2; 
qy(e)=Floadj(ev)-M(ev.ev)*xddj(ey): 
qz(e)=Floadj(ez)-M(ez,ez)*xddj(ez): 
% 
% Vertical or horizontal shear force at element e 
Vy(e)=sum(qy): 
Vz(e)=sum(qz): 
% 
% Vertical or horizontal bending moment at element e 
Myl(e)=Vy(e)*Lece): 
Mzl(e)=Vz(e)*Le(e): 
My(e)=sum(Myl): % BM about z-axis 
Mz(e)=sum(Mz1): % BM about y-axis 
% 
% Axial stress at crown of element e (positive is tension) 
sıgc(e)=-My(e)*Yc(e)/laz(e): 
% Axial stress at keel of element e 
sıgk(e)=My(e)*Yk(e)/laz(e): 
% Axial stress at port-side fiber (port=pos z) (positive is tension) 
sigp(e)=-Mz(e)*Zp(e)/lay(e): 
% Axial stress at stbd-side fiber 
sigs(e)=Mz(e)*Zs(e)/lay(e): 
% 
% Axial strains 
EPSc(e)=sigc(e)/E; 
EPSk(e)=sigk(e)/E: 
EPSp(e)=sigp(e)/E: 
EPSs(e)=sigs(e)/E: 
% 
end 
% Save current time. displacements and velocities 
% (at specified nodes) and axial strains 
% (at specified elements) for later plotting... also save 
% current horizontal fluid velocity and acceleration, and 
% (at nodes 10 and 20) and all calculated honzontal forces 
% (at nodes 10 and 20) for later plotting. 
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ume(1)-t; 
% Vertical/horizontal displacements/velocities 
% at specified nodes 
for n=1:4 
Dv(n)=(N(n)-1)*4+1: 
Dh(n)=(N(n)-1)*4+2: 
end 
xv 1()=xj(Dv(1)): 
xh1(i)=xj(Dh(1)): 
xdv 1(1)=xdj(Dv(1)): 
xdhl(i)=xdj(Dh(1)): 
xv2 (1) -NX( Dv(2)): 
xh2(1)-Nj( Dh(2)): 
xdv2(1)=xdj(Dv(2)): 
xdh2(1)-xdj(Dh(2)); 
xv3(i) -Xj( Dv(3)): 
xh3(1)=xj(Dh(3)): 
xdv3(1)-xdj(Dv(3)): 
xdh3(1)=xdj(Dh(3 )): 
xV4(1)=x}(Dv(4)): 
xh4(i)=sj(Dh(4)): 
xdv4(i)=xdj(Dv(4)): 
xdh4(1)=xdj(Dh(4)): 
% 
xabs(1)=xaj: 
xdabs(1)=xadj: 
xddabs(1)=xadd)j: 
% 
% Horizontal hull/fluid motions and 
% external forces 
xddh2(1)=xddj(Dh(2)): 
xddh3(i)-xddj(Dh(3)); 
uh3(1)=uz(N(3)): 
udh3(1)=uzd(N(3)): 
Fvelh3(1)=Fvel(Dh(3)): 
Facch3(1)=Facc(Dh(3)): 
Floadh3(1)=Floadj(Dh(3)): 
% 
xddh4(1)=xddj(Dh(4)); 
uh4(i)=uz(N(4)); 
udh4(1)=uzd(N(4)): 
Fvelh4(1)=Fvel(Dh(4)); 
Facch4(1)=Facc(Dh(4)); 
Floadh4(1)=Floadj(Dh(4)): 
% 
% Axial strains (keel, crown, port, stbd) at 
% at specified elements 
straink1(1)=EPSk(Ell): 
strainc l(i)>EPSc(EIl): 
strainpl(1)=EPSp(EIl): 
strainsl(1)=EPSs(Ell); 
straink2(1)>=EPSK(El y; 
strainc2(1)=EPSc(El2): 
strainp2(i)=EPSp(EI2): 





end 


Sstrains2(1)=EPSs(E12). 

straink3(1)=EPSk(E13). 

strainc3(1)=EPSc(E]3): 

strainp3(1)- EPSp(E13). 

strains3(1)=EPSs(El3 ); 

% 

% Save the current time step displacements. velocities. and 
% accelerations. and forces for use in the next time step 
XIEN]. 

xdt=xdj: 

xddt-xddj: 

Floadt-Floadj: 

xat-xaj: 

xadt=xadj: 

xaddt-xaddj: 

break. end 


if j==jlimit. break, end 


end 
if j==jlimit, break, end 
if t>=tf. break. end 
end 
% 





% FLUID.M 

% 

% Calculates the free field fluid velocities and accelerations (at each node) 

% 

function[uv.uz.uvd.uzd]-fluidvel(xj.xchg.vchg.zchg.Dship.t. TM.BR.BD.BRD.BDD... 

....BRDD.BDDD.DEL.L1.T1.btime.bdpth.NDOF.NNODES.Lh) 

% 

% Non-dimensional bubble parameters 

Tl=UTI1: 

a=interpl(TM.BR.TI).: 

s=interpl(TM.BRD. TI). 

sd=interp1(TM.BRDD.TI). 

|-interpl (TM.BDD.TIT). 

Id=interpl(TM.BDDD. TI): 

d=interpl(TM.DEL.TD: 

% Source strengths 

e1=(L143)*(a%2)*s/T 1: 

e2=-(L1%4/(2*T1))*(a%3* 1): 

Big IT 1^3/T 1^2)*(a^2*sd-*-2*a*s^2): 

E2d—-(L1^4/(2*T 1^2))*(a^3*ld--3*a^2*1*s); 

% Vertical bubble velocity 

vm=-1*L1/T1; 

% Bubble depth 

ti=t: 

bdepth=interp 1 (btime.bdpth., ti): 

% 

for n=l: NNODES 
ny=(n-1)%4+1: 
nz=(n-1)*4+2; 
xnode=(n-1)*Lh(n); 
ynode-xj(nv); 
znode-xj(nz); 
% vr and hr (vertical and horiz relative distance) and angle of flow 
vr=bdepth-Dshiptynode: 
hr-((xnode-xchg)^2-*(znode-zchg)^2)^0. 5; 
theta=atan((xnode-xchg)/(znode-zchg)): 
% 
% Free-field fluid velocities/accelerations (vertical and horizontal) 
% Horizontal and vertical free-field fluid velocities and accelerations 
dis=(vr42+hr’2)*0.5: 
uhz(el *hr/dis^3)*(3*e2*vr*hr/dis^3). 
uv-(el*vr/dis^3)*(3*e2*vr^2/dis^3)-(e2/dis^3).; 
uhd-(eld*hr/dis^3)*-(3*e2d*hr*vr/dis^5)-... 
...Vm*(3*el *hr*vr/dis^5-(3*e2*hr/dis^5)*(1-(3*vr^2/dis^2))): 
uvd-(eld*vr/dis^3-*3*e2d*vr^2/dis^5-e2d/dis^3)-... 
...vm*(el/dis^3-3*el*vr^2/dis^5-9*e2*vr/dis^5-13*e2*vr^3/dis^7); 
% 
% Y and Z directed free-field fluid velocities and accelerations 
uy(n)-uv; 
uz(n)=uh*costtheta): 
uyd(n)=uvd; 
uzd(n)=uhd*cos(theta); 

end 


p 
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"DRY" WHIPPING RESPONSE AT SPECIFIED NODES 
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Figure B-1. Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode | flexure. 
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Horizontal displacement vs. time for node 10 
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Figure B-2. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode | flexure. 
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Figure B-3. Horizontal response at node 20, given initial displacments equal to a scaled 


horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 7 
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Figure B-4. Horizontal response at node 7, given initial dısplacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Figure B-5. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Figure B-6. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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APPENDIX C 


"WET" WHIPPING RESPONSE AT SPECIFIED NODES 
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Figure C-1. Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 10 
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Figure C-2. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode |} flexure. 
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Horizontal displacement vs. time for node 20 
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Figure C-3. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 1 flexure. 
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Horizontal displacement vs. time for node 7 
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Figure C-4. Horizontal response at node 7, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Horizontal displacement vs. time for node 10 
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Figure C-5. Horizontal response at node 10, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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Figure C-6. Horizontal response at node 20, given initial displacments equal to a scaled 
horizontal mode shape vector for mode 2 flexure. 
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APPENDIX D 


WHIPPING RESPONSE FOR HULL SUBJECTED TO CHARGE DESIGNED TO 
EXCITE HORIZONTAL MODES 1 AND 2 FLEXURAL WHIPPING 
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Figure D-1. Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite horizontal mode 1 flexural whipping. 
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Figure D-2. Horizontal response at node 10, for hull subject to bubble pulse designed to 
excite horizontal mode 1 flexural whipping. 
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Horizontal displacement vs. time for node 20 
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Figure D-3. Horizontal response at node 20, for hull subject to bubble pulse designed to 
excite horizontal mode 1 flexural whipping. 
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Horizontal displacement vs. time for node 7 
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Figure D-4. Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite horizontal mode 2 flexural whipping. 
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Horizontal displacement vs. time for node 10 
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Figure D-5. Horizontal response at node 10, for hull subject to bubble pulse designed to 
excite horizontal mode 2 flexural whipping. 
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Horizontal displacement vs. time for node 20 
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Figure D-6. Horizontal response at node 20, for hull subject to bubble pulse designed to 
excite horizontal mode 2 flexural whipping. 
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APPENDIX E 


WHIPPING RESPONSE FOR HULL SUBJECTED TO CHARGE DESIGNED TO 
EXCITE HORIZONTAL MODE 2 FLEXURAL WHIPPING AND INCLUDING 
INTERNAL MASS/ABSORBER 


53 





z 
ER 


P 
= 


B^. 


Horizontal displacement vs. time for node 7 





0.4 
To? 
(= 
E 
a, O 
O 
o 
2 -0.2 
C 
-0.4 | 
0 0.5 1 55 2 25 3 
Time (sec) 
Horizontal velocity vs. time for node 7 
10 
3 5 
EX 
O 
9 
2 -5 
-10 
0 05 1 155 2 249 3 
Time (sec) 
Horizontal acceleration vs. time for node 7 
— 500 
N 
O 
(D 
m 
€ 
© 0 
© 
® 
® 
O 
< 
-500 
0.5 1 TS 2 25 3 
Time (sec) 


Figure E-1. Horizontal response at node 7, for hull subject to bubble pulse designed to 
excite horizontal mode flexural whipping and including 200 LT "tuned" mass damper at 
node 7. 
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Horizontal displacement vs. time for node 10 
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Figure E-2. Horizontal response at node 10, for hull subject to bubble pulse designed to 
excite horizontal mode flexural whipping and including 200 LT "tuned” mass damper at 
node 7. 
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Horizontal displacement vs. time for node 20 
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F igure E-3. Horizontal response at node 20, for hull subject to bubble pulse designed to 
excite horizontal mode flexural whipping and including 200 LT "tuned" mass damper at 
node 7. 





Horizontal displacement vs. time for internal absorber 
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Figure E-4. Horizontal response of 200 LT "tuned" mass located at node 7. 





APPENDIX F 


NATURAL SLOSHING FREQUENCIES AND SLOSHING MODAL MASS FOR 
LATERAL SLOSHING OF RECTANGULAR TANKS 
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Figure F-1. Natural sloshing frequency and sloshing modal mass for lateral sloshing of a 
(10'h x 10'w x 10'l) rectangular tank. 
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Figure F-2. Natural sloshing frequency and sloshing modal mass for lateral sloshing of a 
(20'h x 10'w x 20'l) rectangular tank. 
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APPENDIX G 


STRAIN AND VELOCITY HISTORIES (MEASURED AND CALCULATED) FOR 
THE "RED SNAPPER" MODEL TESTS 
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Figure G-1. Starboard-side fiber strain histories for elements 6 and 9 for Test 2 of 
the "Red Snapper" model test (horizontal and vertical mode | whipping), and strains 
calculated using the computational algorithm. 
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Figure G-2. Crown and keel fiber strain histories for elements 6 and 9 (respectively) for 
Test 2 of the "Red Snapper" model test (horizontal and vertical mode 1 whipping), and 
strains calculated using the computational algorithm. 


169 





Horizontal velocity vs. time, node 10 
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Figure G-3. Horizontal velocity histories for nodes 10 and 20 for Test 2 of 
the "Red Snapper" model test (horizontal and vertical mode 1 whipping), and velocities 
calculated using the computational algorithm. 
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Axial strain vs. time, element 6 
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Figure G-4. Starboard-side fiber strain histories for elements 6 and 17 for Test 3 of 
the "Red Snapper" model test (horizontal mode 2 whipping), and strains calculated using 
the computational algonthm. 
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Horizontal velocity vs. time, node 10 
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Figure G-5. Horizontal velocity histories for nodes 10 and 20 for Test 3 of 
the "Red Snapper" mode! test (horizontal mode 2 whipping), and velocities calculated 
using the computational algorithm. 
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